t 


t 




AD 


Grant Number DAMD17-94-J-4384 


TITLE: New Methods for Quantitative, High-Resolution Ultrasonic 

Imaging of the Breast 


PRINCIPAL INVESTIGATOR: Robert C. Waag, Ph.D. 

CONTRACTING ORGANIZATION: University of Rochester 

Rochester, New York 14627 


REPORT DATE: January 1998 


TYPE OF REPORT: Final 


PREPARED FOR: U.S. Army Medical Research and Materiel Command 

Fort Detrick, Maryland 21702-5012 


DISTRIBUTION STATEMENT: Approved for public release; 

distribution unlimited 


The views, opinions and/or findings contained in this report are 
those of the author(s) and should not be construed as an official 
Department of the Army position, policy or decision unless so 
designated by other documentation. 


19980617 



§¥t§ INSPECTED 1 


1 




REPORT DOCUMENTATION PAGE 

Forth Approved 

OMB No. 0704-0188 

♦Public reporting burden for this collection of Information Is estimated to average 1 hour per response, including tire time for reviewing instructions, searching existing data sources 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Sendrcomments regarding this burden estimate or any other aspect of this' 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate tor Information Operations and Reports, 121b Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget Paperwork Reduction Project (0704-0188), Washington, DC 20503 

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 

January 1998_ 

3. REPORT TYPE AND DATES COVERED 

Final f30 Sen 94 - 31 Dec 971 

■ - _ ■ ^ ■ - *--*■ — * 1 - 

4. TITLE AND SUBTITLE 

5. FUNDING NUMBERS 
_DAMD17-94-J-4 3R4 

New Methods tor Onanhi hat l ve. Hi nh— Resol nt i on FIT trpsoni r 

Imaging of the Breast 




6. AUTHOR(S) 

Robert C- Waag,.... Ph. D. 





7. PERFORMING ORGANI7ATION NAMRS1 ANn Annncccf^S) 

University of Rochester 

Rochester, New York 14627 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 






9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

U.S. Army Medical Research and Materiel Command 

10. SPONSORING/MONITORING 

AGENCY REPORT NUMBER 

Eoxt .Detrick, Maryland 21702-5012 







11. SUPPLEMENTARY NOTES 





12a. DISTRIBUTION / AVAILABILITY STATEMENT 

12b. DISTRIBUTION CODE 

Approved for public release; distribution unlimited 








13. ABSTRACT (Maximum 200 

The objective of this research is to form quantitative, high-resc 
of the breast using advances in adaptive focussing, a b 
scattering, and unique experimental apparatus. Two new 
implemented and test images were formed to demonstrate the 
to improve resolution for better diagnosis of breast cancer, 
adaptive focussing to compensate for aberration in pulse-ech 
technique, transmit signals are synthesized from scatte 
backpropagation followed by time-shift estimation and conn 
method uses eigenfunctions and eigenvalues of an operator ae 
measurements to focus adaptively on distributed scattering 
quantitative images efficiently. Image evaluation employed co 
morphology of scattering objects as well as with conventioi 
reconstructed x-ray images. 

>lution ultrasonic images 
reakth rough in inverse 
imaging methods were 
potential of the methods 

One method employs 
o imaging. In this new 
red signals by using 
ipensation. The other 
>sociated with scattering 

3 objects and to form 
mparison with observed 
lal b-scan images and 

14. SUBJECT TERMS 

Breast Cancer 

Ultrasound, imaging, diagnosis, quantitative, high-resolution. 

15. NUMBER OF PAGES 

97 : 

16. PRICE CODE 

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 
OF REPORT OF THIS PAGE OF ABSTRACT 

Unclassified- Unclassified Unclassified 

20. LIMITATION OF ABSTRACT 

Unlimited 


NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89) 


Prescribed by ANSI Std. Z39-18 
298-102 


2 





' FOREWORD 

* 


Opinions, interpretations, conclusions and recommendations are 
those of the author and are not necessarily endorsed by the U.S. 
Army. 


_ Where copyrighted material is quoted, permission has been 

obtained to use such material. 

_ Where material from documents designated for limited 

distribution is quoted, permission has been obtained to use the 
material. 

_ Citations of commercial organizations and trade names in 

this report do not constitute an official Department of Army 
endorsement or approval of the products or services of these 
organizations. 

_ In conducting research using animals, the investigator(s) 

adhered to the "Guide for the Care and Use of Laboratory 
Animals," prepared by the Committee on Care and Use of Laboratory 
Animals of the Institute of Laboratory Resources, National 
Research Council (NIH Publication No. 86-23, Revised 1985). 

^ E’or the protection of human subjects, the investigator (s) 
adhered to policies of applicable Federal Law 45 CFR 46. 

_ In conducting research utilizing recombinant DNA technology, 

the investigator(s) adhered to current guidelines promulgated by 
the National Institutes of Health. 

_ In the conduct of research utilizing recombinant DNA, the 

investigator(s) adhered to the NIH Guidelines for Research 
Involving Recombinant DNA Molecules. 

_ In the conduct of research involving hazardous organisms, 

the investigator(s) adhered to the CDC-NIH Guide for Biosafety in 
Microbiological and Biomedical Laboratories. 



3 







4. Table of Contents 

1. Front Cover 1 

2. SF 298 2 

3. Foreword 3 

4. Table of Contents 4 

5. Introduction 5 

6. Body 6 

a. Method 6 

b. Results 6 

c. Discussion 8 

7. Conclusions 9 

8. References 10 

9. Final Report Information 11 

a. Bibliography 11 

b. Personnel 14 

10. List of Appendices 15 


4 



5. Introduction 


> 


Widely recognized strengths of ultrasonic imaging techniques for diagnosis and 
monitoring for recurrence of breast cancer are the nonionizing nature of acoustic 
waves, the ability to provide good contrast between fluid and parenchymal tissue, 
portability, and good patient tolerance. Ultrasonic imaging is inexpensive compared to 
magnetic resonance imaging and is competitive in cost with x-ray mammography. 
However, these strengths are presently offset by significant weaknesses that include 
low resolution, inability to distinguish between solid masses, and the lack of quantitative 
information about breast tissue characteristics. 

For some years, recognition has existed that high-resolution and quantitative ultrasonic 
imaging has the potential to provide diagnostic information not available through 
conventional ultrasonic techniques such as b-scanning. Aberration correction promises 
to increase significantly the resolution of ultrasonic mammography so that small lesions, 
tumors, and microcalcifications can be clearly imaged even through fatty or dense 
breast tissue. Quantitative, high-resolution ultrasonic imaging is particularly desirable 
for the diagnosis of breast disease since the approach has the potential to differentiate 
between types of cancerous tissue, to detect and characterize microcalcifications, to 
detect small tumors even in dense breast tissue, and to yield diagnostically useful 
images without operator-dependent adjustments of parameters. 

The purpose of this research is to develop new ultrasonic imaging concepts based on 
significant advances in adaptive focusing and breakthroughs in inverse scattering as 
well as on newly available novel experimental apparatus. To focus adaptively, transmit 
signals are synthesized from scattered signals by using backpropagation followed by 
time-shift estimation and compensation. To reconstruct quantitatively, eigenfunctions 
and eigenvalues of an operator associated with scattering measurements are used to 
obtain efficiently images of distributed scattering objects. 


5 



6. Body 


a. Method 

The method employed to achieve the goal of quantitative and high-resolution ultrasonic 
imaging employs a combination of hardware innovations, careful measurements of 
ultrasonic aberrations and scattering, and implementation of new ideas for aberration 
correction as well as for image reconstruction via inverse scattering. The novel 
hardware consists of two different systems. One is a custom made pulse-echo 
apparatus in which the size of the transmit aperture can easily be altered and in which 
scattering objects can be easily changed without disturbing the aberrator. The other 
apparatus consists of a ring transducer array with 2048 elements and associated 
circuitry for multiplexing and signal control. Measurements performed with this ring 
produce experimental data from which the degrading effects of wavefront distortion 
imaging were studied. Measurements also provided experimental data from phantoms 
with known acoustic parameters for image reconstruction. Receive beamforming used 
backpropagation and other methods of space-time processing. Measured scattering 
data with three hundred and sixty degrees of spatial coverage was used for the 
reconstruction of quantitative ultrasonic images through the use of newly developed 
inverse scattering methods. 

b. Results 

Ultrasonic wavefront distortion produced by transmission through breast tissue 
specimens was measured in a two-dimensional aperture. Differences in arrival time 
and energy level between the measured waveforms and references that account for 
geometric delay and spreading were calculated as was the waveform similarity factor. 
Details of this investigation were reported in the Journal of the Acoustical Society of 
America. A copy of this paper is included as Appendix A. In this paper, maps of 
distortion are shown in Figure 2 focussing results are shown in Figure 3 and 
comparisons with corresponding abdominal wall data are given in Figures 7 and 8. 

A new method of waveform design was developed and implemented to produce a 
spatially limited plane wave for illumination of the scattering object to be imaged. The 
method employed backpropagation of the desired spatially limited plane wavefront to 
obtain excitation signals for the ring transducer elements. A paper describing the 
methods in detail was published in the IEEE Transactions on Ultrasonics, 
Ferroelectrics, and Frequency Control. A copy of this paper is included as Appendix B. 
Figure 6 in this paper shows illustrative results. 

Estimation and correction of ultrasonic wavefront distortion using pulse-echo data has 
also been implemented in an experimental setting. The processing used random 
scattering measured in a specially designed pulse-echo configuration in which the 
transmit beam is unperturbed by the aberrating medium. Major results in are: 
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1) A tightly focussed transmit beam, such as that achieved with f/1.5, is needed to 
produce a sufficiently coherent scattered wavefront so that time delay can be 
accurately estimated. 

2) Accommodation of waveform shape distortion in time-delay estimation, such as in 
the modified least-mean-square-error method that incorporates a global smoothness 
requirement, is necessary for estimation of the time-delay surface of a wavefront 
that has propagated through distributed inhomogeneities. 

3) The isoplanatic path size usually decreases with increasing severity of distortion and 
is influenced by compensation technique. 

4) Incorporation of waveform shape compensation, such as by using backpropagation 
prior to time-shift compensation, improves the contrast ratio over that obtained by 
time-shift compensation alone. 

A paper that details this work has been accepted for publication in the IEEE 
Transactions on Ultrasonics, Ferroelectrics, and Frequency Control and a preprint of 
this paper is included as Appendix C. Figures 9-14 of the preprint show results that are 
representative. 

A novel inverse scattering method that uses eigenfunctions of the scattering operator 
was examined theoretically and using computations. A unified framework that 
encompasses eigenfunction methods of focussing and image reconstruction in arbitrary 
media was developed. A paper detailing this work has been published in the Journal of 
the Acoustical Society of America and is included as Appendix D to provide more detail. 
Figures 2 through 8 in this paper show representative images. 

Representative high-resolution b-scan images produced by our ring transducer system 
are presented in Figure 1 along with a corresponding b-scan image produced by a 
conventional ultrasonic scanner and a corresponding x-ray image produced by 
computed tomography. The results show that our ring transducer system is capable of 
producing images that are superior to those of commercial b-scan ultrasonic imagers 
and comparable to those obtained from computed x-ray tomography. 
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Fig. 1. Images of an Anthropomorphic 
Breast Phantom, The phantom was 
imaged in a plane containing two 9 mm 
diameter cysts on the periphery of the 
phantom and one 4 mm diameter cyst 
near the center. Top Left: single b-scan 
image obtained with the ring transducer 
system using f/1.0 aperture and scan 
plane diameter of 80 mm. Top right: b- 
scan image produced by a state-of-the-art 
clinical b-scanner operating at the same 
center frequency as the ring transducer 
system. Bottom Left: compounded b- 
scan image from six geometrically 
corrected views with the ring transducer 
system. Bottom Right: X-ray computed 
tomograph obtained with a state-of-the- 
art scanner using a 1 mm slice thickness. 
These images show that ring transducer 
b-scan resolution is superior to current 
clinical scanners at the same frequency 
and is comparable with x-ray computed 
tomography. 


c. Discussion 

The studies summarized above and detailed in appendices show that all of the tasks in 
the statement of work have been largely accomplished although studies employing 
breast tissue and comparative evaluations of imaging modalities were more limited than 
originally planned. These limitations resulted from lost time caused by unforeseen 
equipment failures and the need for more software development than made possible by 
the budget. Nevertheless, the results demonstrate the new ultrasonic imaging methods 
developed in this project have potential. 






7. Conclusions 


The forgoing results and discussion show that the research generally progressed along 
the lines described in the original proposal. Unanticipated problems limited the depth of 
some studies and the extent of progress but important theoretical, experimental, and 
computational results have been obtained. Additional studies that involve more use of 
tissue and further comparative evaluation of imaging modalities are needed to 
determine the ultimate capability of the imaging techniques investigated in this project 
to improve detection, diagnosis, and treatment monitoring of breast cancer. 
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Measurement and correction of ultrasonic pulse distortion 
produced by the human breast 
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Ultrasonic wavefront distortion produced by transmission through breast tissue specimens was 
measured in a two-dimensional aperture. Differences in arrival time and energy level between the 
measured waveforms and references that account for geometric delay and spreading were 
calculated. Also calculated was a waveform similarity factor that is decreased from 1.0 by changes 
in waveform shape. For nine different breast specimens, the arrival time fluctuations had an average 
(±s.d.) rms value of 66.8 (±12.6) ns and an associated correlation length of 4.3 (±1.1) mm, while 
the energy level fluctuations had an average rms value of 5.0 (±0.5) dB and a correlation length of 
3.4 (±0.8) mm. The corresponding waveform similarity factor was 0.910 (±0.023). The effect of 
the wavefront distortion on focusing and the ability of time-shift compensation to remove the 
distortion were evaluated by comparing parameters such as the — 30-dB effective radius, the 
— 10-dB peripheral energy ratio, and the level at which the effective radius departs from an ideal by 
10% for the focus obtained without compensation, with time-shift estimation and compensation in 
the aperture, and with time-shift estimation and compensation performed after backpropagation. For 
the nine specimens, the average —10-dB peripheral energy ratio of the focused beams fell from 3.82 
(±1.83) for the uncompensated data to 0.96 (±0.18) with time-shift compensation in the aperture 
and to 0.63 (±0.07) with time-shift compensation after backpropagation. The average -30-dB 
effective radius and average 10% deviation level were 4.5 (±0.8) mm and —19.2 (±3.5) dB, 
respectively, for compensation in the aperture and 3.2 (±0.7) mm and —22.8 (±2.8) dB, 
respectively, for compensation after backpropagation. The corresponding radius for the 
uncompensated data was not meaningful because the dynamic range of the focus was generally less 
than 30 dB in the elevation direction, while the average 10% deviation level for the uncompensated 
data was -4.9 (±4.1) dB. The results indicate that wavefront distortion produced by breast 
significantly degrades ultrasonic focus in the low MHz frequency range and that much of this 
degradation can be eliminated using wavefront backpropagation and time-shift compensation. 

PACS numbers: 43.80.Cs, 43.80.Ev, 43.80.Vj 


INTRODUCTION 

Widely recognized strengths of ultrasonic imaging tech¬ 
niques for diagnosis and monitoring of breast disease are the 
nonionizing nature of acoustic waves and the ability to pro¬ 
vide good contrast between fluids and parenchymal tissues. 
However, despite advances in transducer technology, breast 
ultrasonography has thus far been relegated to ancillary use, 
largely because resolution is inadequate. 1,2 This has led to 
consideration of the limitations imposed on ultrasonic imag¬ 
ing of the breast by wavefront distortion that arises from 
propagation through breast tissue. 

Several researchers have studied ultrasonic distortion 
produced by the breast and considered how this distortion 
may be compensated. An early in vivo pulse-echo study by 
Moshfeghi and Waag 3 showed that increasing the aperture to 
//1.0 from //2.6 for breast produced only about \ the resolu¬ 
tion improvement predicted in a homogeneous medium. Tra- 
hey et al 4 made one-dimensional transmission measure¬ 


ments of ultrasonic phase distortion caused in vivo by 
propagation through breast and found an average rms arrival 
time aberration of 36.1 ns for 22 subjects. The same group 
later extended their transmission measurements to two di¬ 
mensions and obtained an rms value of 55.3 ns for seven 
volunteers. 5 They concluded that phase distortion should be 
measured and corrected in two dimensions but did not men¬ 
tion amplitude distortion. In other in vivo one-dimensional 
transmission measurements, Zhu and Steinberg observed se¬ 
vere amplitude distortion, which they attributed to refraction 
in addition to scattering, 6,7 and obtained a relation between 
the average sidelobe floor and the normalized variance of the 
amplitude distortion produced by the breast. 8,9 They con¬ 
cluded that large two-dimensional arrays and new algorithms 
that correct both phase and amplitude distortion in two di¬ 
mensions may be needed to reduce the effects of distortion 
produced by the breast. 

Two basic algorithms for correction of ultrasonic distor- 
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tion have been investigated for general imaging applications. 
The first relies on cross correlation of signals in the aperture 
for estimation of pulse arrival time, 10,11 while the other ad- 
justs beamformer delays to maximize signal brightness. ’ 
However, neither of these methods addresses the problems of 
amplitude and waveform distortion, since both use only 
time-shift compensation in the aperture. Another technique 
has been investigated for the removal of amplitude and 
waveform distortion as well as time-shift distortion in spe¬ 
cialized applications when a point source is present as may 
be the case in lithotripsy. This technique employs time- 
reversed signals but is limited because a suitable point source 
is not generally available in every isoplanatic patch to be 
imaged. 14,15 A more recent method developed by Liu and 
Waag 16 can be used to remove time-shift, amplitude, and 
waveform distortion in general imaging applications. Their 
approach models the distorting medium as a phase screen 
placed some distance from the receiving aperture and re¬ 
moves amplitude and waveform distortion by backpropaga- 
tion of the wavefront before applying time-shift compensa¬ 
tion. 

This paper reports a study of ultrasonic wavefront dis¬ 
tortion produced by breast and the effectiveness of the back- 
propagation method in removing the distortion. In the study, 
wavefronts perturbed by transmission through breast tissue 
were measured in a two-dimensional aperture. Statistics de¬ 
scribing arrival time variations, energy level fluctuations, and 
wave shape distortion were calculated and compared to val- 
ues from analogous measurements using abdominal wall. 
The received waveforms were then focused without compen¬ 
sation, with time-shift estimation and compensation in the 
measurement aperture, and with backpropagation followed 
by time-shift estimation and compensation to determine the 
effectiveness of time-shift compensation with and without 
wavefront backpropagation for the improvement of focusing. 

I. METHOD 

Breast tissue specimens were obtained fresh from reduc¬ 
tion mammoplasty surgery and were stored frozen if not used 
immediately for measurements. The specimens came from 
regions of the breast away from the nipple and consisted of 
fat, glandular and connective tissue, and a surface covering 
of skin. Each specimen was essentially planar and had a 
surface area of at least 7X11 cm 2 . The average thickness was 
26.9 mm. The tissue donors were women ranging in age 
from 18 to 65 with a mean age of 34 years. 

Measurements were carried out using the procedure de¬ 
tailed in Ref. 17 and summarized here for convenience. A 
breast tissue specimen was placed in the specimen holder 
with the skin facing the direction of the receiving transducer 
and pressurized to 500 psi for 30 min in order to dissolve any 
gas bubbles present in the tissue. The specimen holder was 
then mounted in the experimental chamber, which was main¬ 
tained at 37 ±1 °C throughout the measurement. Ultrasonic 
pulses emitted by a hemispheric transducer were received by 
a 128-element linear array immediately after propagation 
through the specimen. A two-dimensional area was scanned 
by translating the array in the elevation direction using an 
automated stage. At each elevation, the array elements were 
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accessed sequentially by a multiplexer. The signal from each 
element was digitized into 12-bit samples for a period of 11.8 
/ms at a rate of 20 MHz. The signal was recorded 19 times at 
each element to permit noise reduction through signal aver¬ 
aging. The nominal center frequency was 3.75 MHz for each 
transducer and the -6-dB bandwidth of the received pulse 
was about 2.2 MHz. The element pitch in the receiving trans¬ 
ducer was 0.72 mm and a reflecting mask reduced the receiv¬ 
ing elevation to 1.44 mm. A period of about 35 min was 
required to record the signals from all 128 array elements at 
each of 32 elevations spaced 1.44 mm apart for a total of 
4096 positions over a 92.16X46.08 mm 2 aperture. The total 
source-receiver separation was about 165 mm and the 
specimen-receiver gap was about 8 mm. 

Three groups of measurements were made. The primary 
set consisted of independent measurements using nine differ¬ 
ent specimens. In addition, two specimens were employed 
for sequential measurements in which a 1-cm layer was re¬ 
moved from the bottom of the specimen before each subse¬ 
quent measurement. For two other specimens, a pair of mea¬ 
surements was made with the source in each of two positions 
located 12 mm apart in the array direction. 

Variations in pulse arrival time were calculated using the 
reference waveform method for arrival time estimation and 
differences in shape among the received waveforms were 
quantified using the waveform similarity factor as described 
in Ref. 16. Energy level fluctuations were calculated as de¬ 
scribed in Ref. 17. The calculations are outlined here to iden¬ 
tify the main steps. A reference waveform was constructed 
for a set of data as an average of all the waveforms meeting 
a cross-correlation criterion for similarity. The reference 
pulse was then cross correlated with each of the original 
waveforms and an arrival time surface was calculated from 
the peaks of the correlation functions. Smoothing was done 
to remove questionable outlying points. A two-dimensional, 
fourth-order polynomial fit to the estimated arrival times was 
used to position a window on the original waveforms, and 
arrival time estimation was repeated using the windowed 
data. Geometric effects were removed by fitting a two- 
dimensional, fourth-order polynomial to the newly calculated 
arrival time surface and subtracting the result to obtain the 
arrival time fluctuations. Fluctuations in energy level across 
the received wavefront were calculated by summing the 
squared amplitudes of the windowed signals at each position, 
converting the results to decibel units, and subtracting a fit¬ 
ted two-dimensional, fourth-order polynomial from the re¬ 
sult. This computation, however, did not employ the 20-dB 
restriction on the dynamic range as in Ref. 17, since no ex¬ 
ceptionally low-energy values occurred. The waveform simi¬ 
larity factor, which ranges from an ideal value of 1 to a 
minimum of 0 and is insensitive to absolute amplitude and 
arrival time like a correlation coefficient, was computed for 
all the windowed waveforms throughout the aperture. 

Two forms of compensation for wavefront distortion 
were applied to compare their effect on focusing. In one, 
time-shift estimation and compensation were performed di¬ 
rectly in the receiving aperture. In the other, time-shift esti¬ 
mation and compensation were performed after the wave- 
front had been backpropagated to the point of maximum 
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FIG. 1. Representative waveforms compensated for geometric travel time, (a) H 2 0. (b) BRS 10. (c) BRS 8a. (d) BRS 8b. Each column of panels shows 
waveforms at sequential increments of 2.88 mm in elevation across the first half of the 46.08-mm aperture. At each elevation, the horizontal coordinate is the 
array direction and spans a distance of 92.16 mm in 0.72-mm increments while the vertical coordinate is time and spans an interval of 2.0 /xs in 0.05-/XS 
increments. Signal amplitude is shown linearly on a gray scale with the maximum signal in the two-dimensional aperture for each measurement represented 
by white and the corresponding negative value by black. 


waveform similarity. The backpropagation used the angular 
spectrum method as described in Ref. 16. The reference 
waveform method was employed in both forms for time-shift 
estimation. 

The efficacy of the two correction procedures, as well as 
the effect of the original distortion, were evaluated by focus¬ 
ing each set of data at 180 mm via a Fourier transform 
method as described in Ref. 18. Each focus was described by 
the -10, -20, and —30 dB effective radii, which are half the 
cube root of the product of the corresponding effective 
widths in the array, elevation, and time directions. The effec¬ 
tive radius several dB down from the peak is a useful mea¬ 
sure of point resolution while the effective radius many dB 
down from the peak is a measure of contrast resolution. The 
-10-dB peripheral energy ratio, which is the ratio of energy 
outside an ellipsoid bounded by the -10 dB effective widths 
to the energy inside that ellipsoid, was also used to describe 
each focus quantitatively. Since the energy outside a speci¬ 
fied region around the main peak of the focus can be viewed 
as an integration of sidelobe intensity, the peripheral energy 
ratio is another measure of contrast resolution. Additionally, 
each focus was described by the 10% deviation level, which 
is the level at which the effective radius becomes 10% larger 
than that produced by ideal waveforms obtained for the data 


set by replicating its average time-shift compensated wave¬ 
form throughout the aperture. This level provides a measure 
of the degree to which the central region of actual focus is 
ideal. 

II. RESULTS 

Representative sets of waveforms recorded after trans¬ 
mission through breast tissue specimens and corrected for 
geometric arrival time are shown in Fig. 1 with correspond¬ 
ing representative water path waveforms to illustrate the 
range and combination of distortion levels encountered in 
this study. The water path waveforms in (a) show minimal 
arrival time and energy level fluctuation and nearly ideal 
waveform similarity. The tissue path waveforms in (b) ex¬ 
hibit low arrival time and energy level fluctuation and mod¬ 
erate waveform distortion. The waveforms in (c) have mod¬ 
erate arrival time variation, high energy level fluctuation, and 
low wave shape distortion. The waveforms in (d) show high 
arrival time and waveform distortion but moderate energy 
level fluctuation. 

Arrival time fluctuations and energy level fluctuations 
produced by nine different breast tissue specimens are shown 
in Fig. 2 and statistics of these data as well as the waveform 


1960 J. Acoust. Soc. Am., Vol. 97, No. 3, March 1995 


Hinkelman et a/.: Ultrasonic pulse distortion by human breast 1960 











































































































































































FIG. 2. Breast tissue path arrival time and energy level fluctuations for nine 
different specimens, (a) BRS 7. (b) BRS 8a. (c) BRS 8b. (d) BRS 9a. (e) 
BRS 9b. (f) BRS 10. (g) BRS 11. (h) BRS 13. (i) BRS 17. In the left panel 
of each pair, arrival time difference is shown on a linear scale with a maxi¬ 
mum arrival time fluctuation of +150 ns represented by white and a mini¬ 
mum arrival time fluctuation of -150 ns represented by black. In the right 
panel of each pair, energy level fluctuations are shown on a log scale with a 
maximum positive excursion of +10 dB represented by white and a maxi¬ 
mum negative excursion of —10 dB represented by black. In all panels, the 
horizontal coordinate is the array direction and spans a distance of 92.16 
mm in 0.72-mm increments while the vertical coordinate corresponds to 
position of the array in elevation and spans a distance of 46.08 mm with 
points interpolated from measurements at 1.44-mm intervals to produce data 
at 0.72-mm increments. 
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similarity factor are given in Table I for each set of measured 
waveforms. The arrival time and energy level fluctuation pat¬ 
terns produced by a given specimen are roughly similar and 
the specimen-to-specimen variability is not great. Consistent 
with these observations, the effective arrival time and energy 
level fluctuation correlation lengths are similar for each 
specimen while the average waveform similarity factor has a 
small standard deviation. 

Focal plane amplitude at representative instants of time 
are shown in Fig. 3 for typical measured waveforms that 
have been focused without compensation, with time-shift 
compensation in the aperture, and with time-shift compensa¬ 
tion following backpropagation along with the focus ob¬ 
tained with ideal data. The spread of the focus obtained from 
uncompensated waveforms is large relative to the spread of 
the focus obtained from ideal data. The focus of the wave¬ 
forms that have been time-shift compensated is more con¬ 
centrated than that obtained from uncompensated waveforms 
and further concentration is obtained in the focus of wave¬ 
forms that have been time-shift compensated following back- 
propagation, but the focus of each is not as concentrated as 
the focus obtained from ideal waveforms. 

Effective radius curves of the focus obtained from the 
same typical measured waveforms for uncompensated, time- 
shift compensated, backpropagated and time-shift compen¬ 
sated, and ideal cases are shown in Fig. 4. The effective 
radius of the focus obtained without compensation departs 
from the ideal case by 10% at a level —5.2 dB below the 
peak. Time-shift compensation in the aperture reduces the 
10% deviation level to —21.4 dB while time-shift compen¬ 
sation following backpropagation has a 10% deviation level 
of -25.8 dB and has an effective radius that is appreciably 
narrower at levels 30 to 40 dB below the peak. 

The focus obtained with each of the nine different sets 
of data as well as the improvement in focus obtained using 
time-shift compensation in the aperture and using time-shift 
compensation after backpropagation are described by the pa¬ 
rameters in Table II. The data show that time-shift compen¬ 
sation in the aperture and time-shift compensation following 
backpropagation result in an effective radius at the focus that 
is similar at the —10- and — 20-dB levels but that time-shift 
compensation following backpropagation improves the 
—30-dB effective radius substantially over that obtained with 
time-shift compensation in the receiving aperture. This indi¬ 
cates that time-shift compensation with or without back- 
propagation yields about the same point resolution but that 
time-shift compensation after backpropagation improves the 
contrast resolution more than time-shift compensation in the 
receiving aperture does. The — 10-dB peripheral energy ratio 
obtained with time-shift compensation in the aperture aver¬ 
ages about 25% of that without compensation while the 
— 10-dB peripheral energy ratio obtained with time-shift 
compensation following backpropagation is about 16% of 
that obtained without compensation. These peripheral energy 
ratios are another indication that the sidelobe level, and 
therefore the contrast ratio, is improved more by backpropa¬ 
gation followed by time-shift compensation than by time- 
shift compensation alone. 

The arrival time and energy level fluctuations produced 
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TABLE I. Breast tissue path statistics of wavefront distortion produced by nine different specimens. 


Specimen 

number 

Specimen 

thickness 

(mm) 

Arrival time fluctuations 
rms 99.5% Effective 

value value corr. len. 

(ns) (ns) (mm) 

Energy level fluctuations 
rms 99.5% Effective 

value value corr. len. 

(dB) (dB) (mm) 

Waveform 

simUarity 

factor 

7 

15-25 

63.1 

196.2 

5.82 

4.86 

14.21 

3.20 

0.940 

8a 

30-35 

66.5 

222.0 

4.88 

6.08 

16.35 

4.85 

0.926 

8b 

30 

85.6 

252.0 

5.16 

5.02 

12.99 

3.41 

0.869 

9a 

40 

76.5 

225.7 

4.81 

5.65 

15.10 

4.25 

0.914 

9b 

35 

79.5 

276.6 

4.20 

4.98 

13.09 

3.04 

0.883 

10 

20-25 

59.5 

212.8 

3.38 

4.67 

12.75 

2.87 

0.908 

11 

20-25 

70.7 

252.2 

4.73 

4.77 

12.44 

2.75 

0.903 

13 

15-20 

49.3 

160.8 

3.36 

4.72 

13.44 

3.39 

0.930 

17 

20-25 

50.5 

162.3 

2.40 

4.53 

12.46 

2.46 

0.918 

mean 

26.9 

66.8 

217.8 

4.30 

5.03 

13.65 

3.36 

0.910 

s.d. 

7.7 

12.6 

39.9 

1.07 

0.51 

1.33 

0.76 

0.023 


by different thicknesses of tissue for two specimens are 
shown in Fig. 5 and described statistically in Table III. The 
arrival time fluctuations, energy level fluctuations, and wave¬ 
form distortion increase with specimen thickness for all but 
one measurement. The arrival time and energy level fluctua¬ 
tion patterns for the different thicknesses of the same speci¬ 
men are seen to vary considerably. The focus characteristics 
associated with these patterns are given in Table IV. 

The arrival time and energy level fluctuations in the ap¬ 
erture and after backpropagation are shown in Fig. 6 for two 
specimens for each of two source positions 12 mm apart. The 
features in the fluctuation patterns obtained for each speci¬ 
men are recognizable but shifted in the patterns obtained 
after the change in source position. The features are also 
changed but still identifiable after backpropagation. The 
similarity of the patterns is indicated quantitatively by the 
distortion statistics given in Table V and by the focus char¬ 
acteristics presented in Table VI for compensation using the 
time shifts calculated from the waveforms produced with the 
source at the focal position. The data in Table VI also show 
that compensation using the time shifts calculated from 
waveforms produced with the source 12 mm from the focal 
position is less effective than compensation using the time 
shifts from waveforms produced with the source at the posi¬ 
tion of focus. 

III. DISCUSSION 

The specimens studied here were sections of breast tis¬ 
sue rather than whole breasts. They came from unusually 
large breasts and ranged from 15-40 mm in thickness with 
nearly planar surfaces. Thus while the measurements illus¬ 
trate elements of distortion produced by the breast, they may 
not describe conditions typically encountered in the clinic. 
Nevertheless, the information obtained in this study under 
precisely controlled experimental conditions adds to the 
available data about the ultrasonic wavefront distortion pro¬ 
duced by the breast and should be useful in the development 
of imaging instruments with better resolution for improved 
ultrasonic diagnosis of breast disease. 

Refraction at the specimen-water interfaces could con¬ 
tribute to the distortion measured in this transmission study. 
Prior investigations 16,17 have shown that this contribution is 
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negligible for abdominal wall, for which the surfaces are the 
skin and the peritoneum. The breast tissue specimens, how¬ 
ever, had a cut surface on their lower side, so a special ex¬ 
periment was undertaken to investigate the contribution of 
refraction at this surface to the measured distortion. A thick 
specimen (BRS 9b) was selected and two measurements 
were made over approximately the same area, one with the 
skin side up and the other with the skin side down. Because 
of the 30-mm difference in propagation distance from the cut 
surface to the receiver in these measurements, differences in 
the distortion patterns from these measurements should be 
evident if refraction at the cut surface were an appreciable 
source of the distortion. In this experiment, however, the 
arrival time fluctuation rms values changed by less than 9% 
while the energy level fluctuation rms values and all corre¬ 
lation length values differed by less than 3%. These changes 

i n 

are within the range expected from reproducibility studies 
given the uncertainty in the scan position and the differences 
in propagation path caused by inverting the specimen. The 
arrival time and energy level fluctuation maps from the two 
measurements were also remarkably similar. Therefore, the 
cut specimen surface, which is smoothed by the kapton 
membrane during measurements, is not considered to con¬ 
tribute significantly to the distortion measured for the breast 
specimens. 

The data processing used to determine the wavefront 
distortion in this study employed the same calculated fourth- 
order polynomial references for the determination of arrival 
and energy fluctuation as were described in the report 17 of 
wavefront distortion by abdominal wall, but incorporated 
other improvements and extensions. The polynomial fits, as 
noted in Ref. 17, avoid the need for accurate knowledge 
about tissue dimensions and local acoustic characteristics 
and compensate for small misalignments and low spatial- 
frequency variations. The adaptive reference waveform 
method used to estimate arrival time differences in this paper 
is considered more accurate than the least-mean-square error 
method employed previously because, on average, fewer 
than 20% of the cross correlations computed for each ab¬ 
dominal wall data set using the adaptive reference waveform 
method had a peak value below 0.8, while about 30% of the 
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FIG. 3. Focal plane time histories of representative measured data (BRS 7). Each panel shows the bipolar distribution of signal amplitude as a shade of gray 
on a 50-dB log scale for each polarity in the focal plane at an instant of time. In all the panels, the horizontal coordinate is elevation and spans 37.504 mm 
while the vertical coordinate is in the array direction and spans 56.256 mm. The number beneath each panel identifies the (zero-origin) instant of time in the 
128-point interval employed in the temporal Fourier transform. First row: Uncompensated data. Second row: After time-shift estimation and compensation in 
the aperture. Third row: After backpropagation of 40 mm followed by time-shift estimation and compensation. Fourth row: Ideal data. 


peaks were below 0.8 for the least-mean-square error 
method. Also, the reference waveform method gives consis¬ 
tently lower estimates of distortion for water paths and yields 
superior focus characteristics when employed for time-shift 
compensation. Therefore, the reference waveform method 
was used in this study and the data in Ref. 16 were repro¬ 
cessed using the reference waveform method to permit a 
comparison of the breast results discussed here with the ab¬ 
dominal wall data. 

In the nine independent measurements reported here, an 
average (±s.d.) of 2807 (±228) waveforms or 68.5% of the 
4096 received waveforms were similar enough to be incor¬ 
porated into the final reference waveform. The average of the 
maximum value for correlations between each of the win¬ 
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dowed waveforms and the final reference waveform was 
0.868 (±0.018). This figure, calculated using all of the wave¬ 
forms in the aperture before any smoothing, is similar to the 
average cross-correlation value of 0.86 calculated for neigh¬ 
boring waveforms in the 35 data sets selected for analysis by 
Freiburger et al 5 in their two-dimensional study of the 
breast. An average of 758 (±183) or 18.5% of the calculated 
delays were deemed unacceptable because they deviated sig¬ 
nificantly from the overall delay surface. Of these, 495 
(about two thirds) were replaced by delays from new peak 
searches while the remaining 263 (about one third) were cor¬ 
rected by smoothing. 

The representative waveforms in Fig. 1 illustrate the 
variability of the wavefronts from which the characteristics 
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FIG. 4. Effective radius of the focus obtained with representative measured 
waveform data (BRS 7). Unc=uncompensated data. TSC=time-shift com¬ 
pensation in the receiving aperture. BP & TSC=backpropagation followed 
by time-shift compensation. Ideal=replication of a single average waveform 

throughout the receiving aperture. 

■ 

of wavefront distortion have been found. The uniformity of 
the water path waveforms indicates that the error introduced 
by the apparatus is small compared to the fluctuations pro¬ 
duced by the tissue path. The presence of tissue creates 
waveform shape changes which decorrelate the waveforms 
and introduce uncertainty in the estimation of arrival time. 
Improvement of arrival time estimation accuracy, and thus 
the correction of wavefront distortion, necessitates the re¬ 
moval of wave shape distortion and has motivated the intro¬ 
duction of a backpropagation step before time-shift compen¬ 
sation. 

The arrival time and energy level fluctuation maps in 
Fig. 2 for the nine independent breast tissue measurements 
show the range of ultrasonic wavefront distortion patterns 
produced by human breast tissue in this study. These patterns 
have different characteristics from those reported for the 
abdominal wall. For example, the breast arrival time fluctua¬ 
tion maps contain smaller features than those in the abdomi¬ 
nal wall arrival time plots. The backgrounds of the breast 
arrival time and energy level fluctuation maps are also com¬ 
prised of larger, more irregular patches than those in the 
abdominal wall plots. However, as is the case with the ab¬ 
dominal wall patterns, some breast arrival time fluctuation 



FIG. 5. Breast tissue path arrival time and energy level fluctuations for 
reductions of specimen thickness, (a)-(c) BRS 9a. (d) and (e) BRS 11. For 
each specimen, the results are presented from top to bottom in order of 
decreasing thickness. The format and scales in each pair of panels are the 
same as in Fig. 2. 

map features appear to correlate with energy level fluctua¬ 
tions while others do not. 

The bar charts in Fig. 7 compare the arrival time, energy 
level, and waveform distortion statistics for the nine breast 
specimen measurements to those in Ref. 17 for abdominal 
wall after recalculation using the reference waveform 
method. In general, breast tissue appears to cause more dis¬ 
tortion than the abdominal wall. The rms arrival time fluc¬ 
tuations produced by the breast specimens in this study have 


TABLE II. Breast tissue path focus characteristics for nine different specimens. r 6 =effective radius. PER=peripheral energy ratio. Unc=uncompensated 
TSC=time-shift compensation. BP=backpropagation followed by TSC. BP dist—distance of backpropagation for maximum waveform similarity. 


-10 dB r € —20 dB r e -30 dB r e -10 dB PER 10% dev. lev. 


Specimen 

number 

BP dist. 
(mm) 

Unc 

(mm) 

TSC 

(mm) 

BP 

(mm) 

Unc 

(mm) 

TSC 

(mm) 

BP 

(mm) 

Unc 

(mm) 

TSC 

(mm) 

BP 

(mm) 

Unc 

TSC 

BP 

Unc 

(dB) 

TSC 

(dB) 

BP 

(dB) 

7 


1.38 

1.04 

wrm 

3.94 

1.52 

1.49 

♦ • ♦ 

3.77 

2.60 

3.350 

0.766 

0.524 

-5.2 

-21.4 

-25.8 

8a 


1.48 

1.09 


4.39 

1.93 

1.69 


4.66 

3.44 

2.733 

0.905 

0.629 

-1.0 

-12.8 

-17.0 

8b 


1.67 

1.07 

siiiK 

7.60 

1.92 

1.54 


5.41 

3.24 

8.031 

1.280 

0.595 

-0.9 

-17.7 

-22.0 

9a 


2.16 

1.06 

MMM 

5.02 

2.06 

1.55 


4.48 

3.76 

3.966 

1.128 

0.654 

-4.9 

-15.2 

-20.9 

9b 


1.52 

1.12 

Mffil 

5.73 

1.63 

1.57 


6.15 

4.49 

4.943 

1.044 

0.742 

-1.0 

-20.2 

-23.7 

10 


1.17 

1.04 

fr^ 

3.39 

1.59 

1.50 


3.64 

2.54 

3.957 

0.968 

0.592 

-8.0 

-19.4 

-22.2 

11 


2.78 

1.04 




1.47 


4.34 


2.618 

0.962 

0.700 

-3.4 

-22.2 

-24.4 

13 


1.16 

1.05 

H—1 

2.71 

1.54 

1.48 

5.92 

3.74 

2.89 

1.766 

0.712 

0.554 

-6.6 

-20.1 

-23.5 

17 



1.01 


2.83 

1.43 

1.40 

7.52 

4.25 

2.57 

3.038 

0.877 

0.631 

-13.3 

-24.0 

-26.1 

mean 



1.06 


4.59 

1.68 

1.52 

» » » 

4.49 

3.17 

3.822 

0.960 

0.625 

-4.9 

-19.2 

-22.8 

s.d. 





1.59 



• • • 



1.827 

0.176 

0.068 

4.1 

3.5 

2.78 
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TABLE III. Breast tissue path statistics of wavefront distortion for reductions of specimen thickness. 


Arrival time fluctuations Energy level fluctuations 


Specimen 

number 

Specimen 

thickness 

(mm) 

rms 

value 

(ns) 

99.5% 

value 

(ns) 

Effective 

corr. len. 
(mm) 

rms 

value 

(dB) 

99.5% 

value 

(dB) 

Effective 

corr. len. 
(mm) 

Waveform 

similarity 

factor 

9a 

40 

76.5 

225.7 

4.81 

5.65 

15.10 

4.25 

0.914 


30 

58.2 

183.6 

4.86 

5.44 

14.79 

4.06 

0.947 


18 

54.6 

157.6 

5.12 

5.01 

15.31 

2.92 

0.965 

11 

20-25 

70.7 

252.2 

4.73 

4.77 

12.44 

2.75 

0.903 


15 

56.2 

168.1 

4.10 

4.92 

14.67 

2.80 

0.933 


a mean of 66.8 ns, which is about 20% higher than the mean 
of 55.5 ns for the 14 abdominal wall samples. The rms en¬ 
ergy level fluctuations produced by breast in this study have 
a mean value of 5.03 dB, which is about 39% higher than the 
mean of 3.62 dB for the abdominal wall. Wave shape distor¬ 
tion is also greater for the breast studies, which have a mean 
waveform similarity factor of 0.910 versus 0.938 for the ab¬ 
dominal wall. 

The difference in the average thickness of breast and 
abdominal wall specimens may contribute to the difference 
in the average arrival time fluctuation rms values for the two 
types of specimens but appears to be a secondary source of 
the variation in energy level fluctuations. The mean thickness 
of the breast specimens is 25% larger than that of the ab¬ 
dominal walls and the average arrival time fluctuation of the 
breast specimens is about 20% larger than that of the ab¬ 
dominal walls. In addition, the arrival time fluctuation rms 


values measured for specimens of both types having similar 
thicknesses were approximately the same. Thus the greater 
thickness of the breast specimens likely contributed impor¬ 
tantly to the greater arrival time fluctuations observed in this 
study. However, a comparison of similarly thick breast and 
abdominal wall specimens revealed that the energy level 
fluctuation rms values of the breast specimens were larger in 
every case. This suggests that breast tissue produces more 
energy fluctuation per unit thickness than does abdominal 
wall tissue, so that the mean rms energy level fluctuation 
produced by the breast specimens would be greater than that 
of the abdominal walls even without the difference in their 
thicknesses. 

The spatial distributions of arrival time and energy level 
fluctuations caused by breast and abdominal wall tissue dif¬ 
fer as well. For the breast, the effective correlation length of 
arrival time fluctuation has a mean of 4.31 mm, which is 



FIG. 6. Arrival time and energy level fluctuations before and after backpropagation for two source positions 12 mm apart, (a) BRS 8a. (b) BRS 17. For each 
specimen, the four-panel set on the left shows arrival time fluctuations and the four-panel set on the right shows the corresponding energy level fluctuations. 
In each four-panel set, the upper row shows the fluctuations with the source to the right and the lower row shows the fluctuations with the source to the left, 
while the left column shows the fluctuations in the measurement aperture and the right column shows the corresponding fluctuations after backpropagation to 
the distance of maximum waveform similarity. The scales in the panels are the same as in Fig. 2. 
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TABLE IV. Breast tissue path focus characteristics for reductions of specimen thickness. r e =effective radius. PER=peripheral energy ratio. Unc 
=uncompensated. TSC=time-shift compensation. BP=backpropagation followed by TSC. BP Dist=distance of backpropagation for maximum waveform 
similarity. 


— 10 dB r € -20 dB r e -30 dB r e -10 dB PER 10% dev. lev. 


Specimen 

number 

BP dist. 
(mm) 

Unc 

(mm) 

TSC 

(mm) 

BP 

(mm) 

Unc 

(mm) 

TSC 

(mm) 

BP 

(mm) 

Unc 

(mm) 

TSC 

(mm) 

BP 

(mm) 

Unc 

TSC 

BP 

Unc 

(dB) 

TSC 

(dB) 

BP 

(dB) 

9a 


2.16 

1.06 

1.04 



1.55 

... 

4.48 

3.76 

3.966 

1.128 

0.654 

-4.9 

-15.2 

-20.9 



1.14 

1.04 

1.03 

3.66 

1.73 

1.56 

7.20 

3.82 

2.89 

2.372 

0.830 

0.540 

-1.0 

-17.0 

-19.8 



1.56 

1.03 

1.03 


1.49 

1.46 

8.06 

4.21 

2.94 

2.954 

0.673 

0.508 

-1.0 

-22.1 

-25.0 

11 


2.78 

1.04 

1.02 



1.47 

* * * 

4.34 

3.03 

2.618 

0.962 

0.700 

-3.4 

-22.2 

-24.4 



1.36 

1.01 

1.01 

3.88 

1.43 

1.41 

• • • 

2.99 

2.52 

2.779 

0.831 

0.565 

-5.6 

-22.2 

- 24.2 


26% smaller than the 5.78 mm mean calculated for the ab¬ 
dominal wall. The effective correlation length of energy level 
fluctuation for the breast has a mean of 3.36 mm, which is 
45%, larger than the abdominal wall correlation length of 
2.31 mm. Consequently, the ratio of the average arrival time 
fluctuation effective correlation length to the average energy 
level fluctuation correlation length is 2.5 for the abdominal 
wall but only 1.3 for the breast specimens. This shows that 
the arrival time and energy level fluctuation correlation 
lengths for the breast are more similar to each other than are 
those for the abdominal wall. 

The focus obtained with uncompensated breast tissue 
waveforms was predictably worse than that obtained with 
abdominal wall waveforms as is illustrated by the bar charts 
in Fig. 8 for uncompensated data. For example, the mean 
— 10-dB peripheral energy ratio of the focused breast data 
was 3.8, nearly twice the average value of 2.1 computed for 
the abdominal data. Also, although not shown in Fig. 8, the 
mean — 20-dB effective radius for the focused breast data 
was 4.6 mm while the corresponding value was 3.6 mm for 
the abdominal wall data. However, both the average breast 
and abdominal wall data sets first deviate 10% from their 
ideal effective radius curves at about —5 dB. 

Time-shift compensation after backpropagation im¬ 
proves the focus of breast data more than time-shift compen¬ 
sation in the aperture does. However, in neither case is the 
average focus better than that of abdominal wall data after 
similar processing. The bar charts in Fig. 8 illustrate these 
trends. For example, the average 10% deviation level for 
breast data is reduced to -18.4 dB by time-shift compensa¬ 
tion in the aperture and to -22.3 dB by time-shift compen¬ 
sation after backpropagation, but the corresponding values 
for the abdominal wall data are —23.5 and —26.3 dB, re¬ 
spectively. Similarly, time-shift compensation in the aperture 
reduces the average -10-dB peripheral energy ratio to 1.0 


for the breast measurements and time-shift compensation af¬ 
ter backpropagation further reduces it to 0.60, while the cor¬ 
responding values for the abdominal wall measurements are 
0.63 and 0.47. Although both compensation techniques im¬ 
prove the focus of ultrasonic waveforms that have propa¬ 
gated through breast tissue, time-shift compensation is more 
effective when performed after backpropagation. 

The distance of backpropagation for maximum wave¬ 
form similarity, i.e., the position of the equivalent phase 
screen in the propagation model employed here, is indicated 
by the data in Tables I and II (as well as in Tables III—VI) to 
be usually greater than the specimen thickness. However, as 
already noted, a distance of about 8 mm is present between 
the top surface of the specimen and the receiving aperture in 
the measurements. Adding 8 mm to the mean (±s.d.) speci¬ 
men thickness given in Table II yields 34.9 (±7.7) mm while 
the mean backpropagation distance is 36.7 (±8.7) mm. Thus 
the data indicate that the optimum backpropagation distance 
is approximately the sum of the specimen thickness and the 
specimen-receiver separation, as in the case of the abdominal 
wall measurements. 16 The physical origin of this circum¬ 
stance is currently unknown, although observations and ex¬ 
periments noted earlier indicate that the origin is not refrac¬ 
tion at the specimen boundary. Additional studies are needed 
to provide information about the relation between tissue mor¬ 
phology and ultrasonic aberration. 

The plots of arrival time and energy level fluctuations in 
Fig. 5 and the statistics of Table III show that reducing the 
tissue path length decreases the severity of the distortion pro¬ 
duced in transmitted ultrasonic pulses. For both specimens 
9a and 11, the magnitude of the measured arrival time and 
energy level fluctuations decreases as the tissue thickness is 
reduced while waveform similarity generally increases. The 
lack of similarity in the patterns of fluctuations as each speci¬ 
men decreases in thickness is attributed to configurational 


TABLE V. Breast tissue path statistics of wavefront distortion for two source positions 12 mm apart. 


Arrival time fluctuations Energy level fluctuations 


Specimen 

number 

Source 

position 

Specimen 

thickness 

(mm) 

rms 

value 

(ns) 

99.5% 

value 

(ns) 

Effective 

corr. len. 
(mm) 

rms 

value 

(dB) 

99.5% 

value 

(dB) 

Effective 
corr. len. 
(mm) 

Waveform 

similarity 

factor 

8a 

Left 

30-35 

66.5 

222.0 

4.88 

6.08 

16.35 

4.85 

0.926 


Right 

30-35 

65.2 

216.0 

4.17 

6.01 

15.11 

4.79 

0.916 

17 

Left 

20-25 

50.5 

162.3 

2.40 

4.53 

12.46 

2.46 

0.918 


Right 

20-25 

47.0 

150.4 

2.27 

4.36 

11.93 

2.44 

0.909 
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TABLE VI. Breast tissue path focus characteristics for two source positions. The numbers in parentheses are values obtained at one source position with 
waveforms compensated using time-shifts calculated with waveforms from the other source position. r () =effective radius. PER ^peripheral energy ratio. 
Unc=uncompensated. TSC=time-shift compensation. BP=backpropagation followed by TSC. BP Dist=distance of backpropagation for maximum waveform 
similarity. 


Specimen 

number 

Source 

position 

BP dist. 
(mm) 

Unc 

(mm) 

-10 dB r 

TSC 

(mm) 

e 

BP 

(mm) 

Unc 

(mm) 

-20 dB 

TSC 

(mm) 

r e 

BP 

(mm) 

Unc 

(mm) 

-30 dB 

TSC 

(mm) 

r e 

BP 

(mm) 

Unc 

10 dB PER 

TSC BP 

Unc 

(dB) 

10% dev. lev. 

TSC BP 

(dB) (dB) 

8a 

Left 

40 

1.48 

1.09 

(1.13) 

1.06 

(1.12) 

4.39 

1.93 

(2.41) 

1.69 

(2.17) 

• - • 

4.66 

(6.76) 

3.44 

(6.24) 

2.73 

0.91 

(1.75) 

0.63 

(1.35) 

-1.0 

-12.8 

(-7.7) 

-17.0 

(-8.8) 


Right 

40 

1.43 

1.05 

(1.09) 

1.05 

(1.09) 

4.61 

1.82 

(2.47) 

1.55 

(1.87) 

... 

4.86 

(7.20) 

3.74 

(6.23) 

3.17 

1.05 

(1.91) 

0.67 

(1.34) 

-1.0 

-16.4 

(-11.5) 

-20.6 

(-12.0) 

17 

Left 

30 

1.04 

1.01 

(0.99) 

1.00 

(1.00) 

2.83 

1.43 

(1.76) 

1.40 

(1.67) 

7.52 

4.25 

(6.98) 

2.57 

(5.46) 

3.04 

0.88 

(2.29) 

0.63 

(1.75) 

-13.3 

-24.0 

(-17.7) 

-26.1 

(-19.1) 


Right 

30 

1.06 

1.02 

(1.01) 

1.01 

(1.02) 

2.86 

1.53 

(1.75) 

1.46 

(1.64) 

7.07 

4.67 

(...) 

2.31 

(5.92) 

2.79 

0.88 

(2.32) 

0.63 

(1.78) 

-12.3 

-22.2 

(-16.9) 

-28.1 

(-18.4) 


changes that result from removal of tissue slabs, differences 
in position of the specimen, and deformation from pressure 
applied when the specimen was sliced. 

The nine separate breast measurements show a similar 
dependence of distortion on specimen thickness. Linear re¬ 
gression establishes significant relationships between breast 
tissue thickness and both rms arrival time and energy level 
fluctuation values, with correlation coefficients of 0.729 and 
0.724, respectively. The energy level fluctuation correlation 
lengths are somewhat correlated to thickness as well. For the 
abdominal wall specimens, the arrival time fluctuation rms 
and correlation length values are moderately correlated to 
thickness and the waveform similarity factor is closely cor¬ 
related, but energy fluctuation rms values and correlation 
lengths are not correlated to thickness. Since breast tissue 
consists of a heterogeneous arrangement of fat, glandular, 
and connective tissue while the abdominal wall is comprised 
of distinct fat and muscle layers, morphological differences 
are thought to be the reason for the observed differences in 
fluctuations produced by breast and abdominal wall. 

The arrival time fluctuation maps in Fig. 6 and focus 
parameters in Table VI for measured wavefronts with differ¬ 
ent source positions demonstrate that the time shifts that best 
compensate the focus at one position may be less effective 
when the desired focus is moved 12 mm laterally and indi¬ 
cate that the use of backpropagation reduces this problem. 
Compensation in the aperture using the time-shifts calculated 
for the adjacent source position, i.e., cross compensation, 
improves the focus in every case, but is far less effective than 
compensation using the time-shifts calculated for the data set 
being corrected. Backpropagation increases the peak cross¬ 
correlation coefficient of the arrival time fluctuation map 
pairs to 0.728 from 0.564 for specimen 8a and to 0.568 from 
0.426 for specimen 17. Cross compensation is more effective 
after backpropagation, but still does not achieve the level of 
focus improvement obtained using regular time-shift com¬ 
pensation in the aperture. This implies a limitation on the 
size of the region over which a single pattern of compensa¬ 
tion is useful. For a related theoretical treatment that yields 
an expression for the size of the isoplanatic patch, see Ref. 
19. 
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The statistics of the arrival time distortion measured in 
this study can be compared to the in vivo results of 
Freiburger et al for a two-dimensional array. 5 Their mean 
rms phase distortion value of 55.3±14.0 ns obtained for 35 
scans on the left breasts of seven volunteers is somewhat 
lower than the 66.8±12.6 ns determined here for nine breast 
specimens. However, their tissue paths had a mean length of 
78.8 ±21.4 mm, 20 which is much longer than those measured 
here. As discussed above, a longer path is expected to pro¬ 
duce larger fluctuations in arrival time. Also, their receiver 
elements were 0.51X3.50 mm 2 , which corresponds to a sur¬ 
face area 1.7 times larger than that of the elements used in 
this study. A larger measurement spot is associated with a 
decrease in measured arrival time fluctuations because the 
fluctuations are averaged over the measurement spot, as de¬ 
tailed in a study 21 that emulated the outputs of elements of 
various sizes in one-dimensional apertures from the mea¬ 
sured outputs of smaller elements in a two-dimensional ap¬ 
erture. Nevertheless, the 4.31±1.08-mm phase correlation 
length measured here is similar to the 4.21±1.14-mm value 
obtained by Freiburger et al. 5 in the azimuthal direction. The 
effects of the distortion measured in the two cases on the 
resulting focus are not comparable because of substantial dif¬ 
ferences in focusing parameters. 

The wavefront distortion and focus degradation found in 
this study may also be compared with the already cited in 
vivo results of Zhu and Steinberg. 6-9 Their observation of 
severe amplitude distortion accompanying phase distortion 
in narrow-band measurements is analogous to energy level 
fluctuations accompanying arrival time fluctuations in the 
present wideband measurements. The emphasis that Zhu and 
Steinberg give to the influence of amplitude distortion in 
their analysis is strongly supported by the new data given 
here. Zhu and Steinberg also showed representative recon¬ 
structions of source profiles that are analogous in the present 
study to effective width curves of the uncompensated focus 
in the array direction. The high degradation of the source 
profiles illustrated in their reports is qualitatively similar to 
the poor focus characteristics given here for uncompensated 
data. However, a quantitative comparison of source profile 
and focus width characteristics is not readily made with the 
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FIG. 7. Comparison of breast and abdominal wall wavefront distortion sta¬ 
tistics. In each chart, the average and standard deviation of the measure¬ 
ments within each group are shown. 


FIG. 8. Comparison of breast and abdominal wall focus statistics. In each 
chart, the average and standard deviation of the measurements within each 
group are shown. Uncomp=uncompensated waveforms. TSC=time-shift 
compensation. BP&TSC=backpropagation before time-shift compensation. 
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published Zhu and Steinberg data because they employed 
different focus descriptors, as well as relatively narrow-band 
pulses, a larger measurement spot size, and significantly 
longer tissue paths. 

IV. CONCLUSION 

Ultrasonic wavefront distortion produced by transmis¬ 
sion through breast specimens has been investigated by re¬ 
cording and processing pulsed waveforms to obtain maps 
and statistics of arrival time and energy level fluctuations. 
The average rms value of the arrival time fluctuations was 
similar to that for abdominal wall specimens of like thick¬ 
ness while the average rms value of energy level fluctuations 
was greater than that for similarly thick abdominal wall 
specimens. The average effective correlation length of the 
arrival time fluctuations for breast specimens was shorter 
than that for the abdominal wall while the average effective 
correlation length of the energy level fluctuations was larger 
than that for the abdominal wall. The breast tissue wave¬ 
forms showed more variation than the abdominal wall wave¬ 
forms. Time-shift compensation in the aperture improved the 
focus of pulses distorted by breast tissue specimens and 
time-shift compensation after backpropagation brought the 
focus closer to the ideal, but the improvement was generally 
less than that obtained for distortion produced by abdominal 
wall. Repeated study of specimens with reduction of thick¬ 
ness showed that pulse arrival times, energy levels, and wave 
shape are increasingly altered as tissue path length increases. 
Measurements using different source positions indicated that 
time-shift compensation in the aperture with arrival times 
estimated when the source was 12 mm away from the focal 
position produced little improvement in the focus, but that 
backpropagation made such cross compensation more effec¬ 
tive, although still less effective than self-compensation in 
the receiving aperture. These results demonstrate that ultra¬ 
sonic propagation through breast tissue produces appreciable 
arrival time and energy level fluctuations and provide impor¬ 
tant new information about the variety and range of these 
fluctuations. They also show that time-shift compensation 
improves the focus of waveforms distorted by breast tissue, 
especially when the compensation is applied after wavefront 
backpropagation. 

ACKNOWLEDGMENTS 

The authors thank Dr. Elethea Caldwell, Dr. Leon Met- 
lay, and Mick Kazee for their assistance in obtaining and 
characterizing breast tissue specimens from surgery. Coop¬ 
eration and assistance from the staff of the Diagnostic Ultra¬ 
sound Department at the University of Rochester’s Strong 
Memorial Hospital is acknowledged with appreciation as are 
many useful discussions with Dr. T. Douglas Mast. Funding 
for this investigation was provided by the University of 
Rochester Diagnostic Ultrasound Research Laboratory In¬ 
dustrial Associates, NIH grants DK 45533 and HL 50855, 


NSF grant BCS92-09680, and U.S. Army Medical Research 

6 Development Command grants DAMD 17-93-J-3014 and 
DAMD 17-94-J-4384. Computations were performed at the 
Cornell National Supercomputing Facility, which is sup¬ 
ported in part by the National Science Foundation, New York 
State, and the IBM Corporation. 

J P M. Jokich, D. L. Monticciolo, and Y. T. Adler, “Breast Ultrasonogra¬ 
phy,” Radiol. Clin. North Am. 30(5), 993-1009 (1992). 

2 W. J. Donegan, “Evaluation of a Palpable Breast Mass,” N. Engl. J. Med. 

327(13), 937-942 (1992). 

3 M. Moshfeghi and R. C. Waag, “ In-Vivo and In-Vitro Ultrasound Beam 
Distortion Measurements of a Large Aperture and a Conventional Aper¬ 
ture Focussed Transducer,” Ultrasound Med. Biol. 14(5), 417-430 
(1988). 

4 G. E. Trahey, P. D. Freiburger, L. F. Nock, and D. C. Sullivan, “/« Vivo 
Measurements of Ultrasonic Beam Distortion in the Breast,” Ultrason. 

Imag. 13(1), 71-90 (1991). 

5 P. D. Freiburger, D. C. Sullivan, B. H. LeBlanc, S. W. Smith, and G. E. 

Trahey, “Two Dimensional Ultrasonic Beam Distortion in the Breast: In 
Vivo Measurements and Effects,” Ultrason. Imag. 14(4), 398-414 (1992). 

6 Q. Zhu and B. D. Steinberg, “Large-Transducer Measurements of Wave- 
front Distortion in the Female Breast,” Ultrason. Imag. 14(3), 276-299 
(1992). 

7 Q. Zhu and B. D. Steinberg, “Wavefront Amplitude Distribution in the 
Female Breast,” J. Acoust. Soc. Am. 96, 1-9 (1994). 

8 Q. Zhu and B. D. Steinberg, “Wavefront Amplitude Distortion and Image 
Sidelobe Levels—Part I: Theory and Computer Simulations,” IEEE Trans. 

Ultrason. Ferroelect. Freq. Control 40(6), 747-753 (1993). 

9 Q. Zhu, B. D. Steinberg, and Ronald Arenson, “Wavefront Amplitude 
Distortion and Image Sidelobe Levels—Part II: In Vivo Experiments,” 

IEEE Trans. Ultrason. Ferroelect. Freq. Control 40(6), 754-762 (1993). 

10 S. W. Flax and M. O’Donnell, “Phase-Aberration Correction Using Sig¬ 
nals from Point Reflectors and Diffuse Scatterers: Basic Principles,” IEEE 
Trans. Ultrason. Ferroelect. Freq. Control 35(6), 758—767 (1988). 

11 M. O’Donnell and S. W. Flax, “Phase-Aberration Correction Using Sig¬ 
nals from Point Reflectors and Diffuse Scatterers: Measurements,” IEEE 
Trans. Ultrason. Ferroelect. Freq. Control 35(6), 768-774 (1988). 

12 L. Nock, G. E. Trahey, and S. W. Smith, “Phase Aberration Correction in 
Medical Ultrasound Using Speckle Brightness as a Quality Factor,” J. 

Acoust. Soc. Am. 85, 1819-1833 (1989). 

13 D. Zhao and G. E. Trahey, “A Statistical Analysis of Phase Aberration 
Correction Using Image Quality Factors in Coherent Imaging Systems,” 

IEEE Trans. Med. Imag. 11(3), 446-452 (1992). ; 

14 M. Fink, “Time Reversal of Ultrasonic Fields—Part I: Basic Principles,” 

IEEE Trans. Ultrason. Ferroelect. Freq. Control 39(5), 555-566 (1992). 

15 F. W. Wu, J.-L. Thomas, and M. Fink, “Time Reversal of Ultrasonic 
Fields—Part II: Experimental Results,” IEEE Trans. Ultrason. Ferroelect. 

Freq. Control 39(5), 567-578 (1992). 

16 D.-L. Liu and R. C. Waag, “Correction of Ultrasonic Wavefront Distortion 
Using Backpropagation and a Reference Waveform Method for Time- 
Shift Compensation,” J. Acoust. Soc. Am. 96, 649-660 (1994). , 

17 L. M. Hinkelman, D.-L. Liu, L. A. Metlay, and R. C. Waag, “Measure- j 

ments of Ultrasonic Pulse Arrival Time and Energy Level Variations Pro- ( 

duced by Propagation through Abdominal Wall,” J. Acoust. Soc. Am. 95, 

530-541 (1994). { 

18 D.-L. Liu and R. C. Waag, “Time-Shift Compensation of Ultrasonic Pulse 
Focus Degradation Using Least-Mean-Square Error Estimates of Arrival 
Time,” J. Acoust. Soc. Am. 95, 542-555 (1994). 

19 B. D. Steinberg and A. K. Luthra, “Simple Theory of the Effects of Me¬ 
dium Turbulence upon Scanning with an Adaptive Phased Array,” J. 

Acoust. Soc. Am. 71, 630-634 (1982). 

20 P. D. Freiburger (personal communication, 20 April 1994). 

21 D.-L. Liu and R. C. Waag, “A Comparison of Wavefront Distortion and 
Compensation in One-Dimensional and Two-Dimensional Apertures,” 

IEEE Trans. Ultrason. Ferroelect. Freq. Control (in press). 


1969 J. Acoust. Soc. Am., Vol. 97, No. 3, March 1995 


Hinkelman et a/.: Ultrasonic pulse distortion by human breast 1969 




Final Report 


Appendix B 


DAMD17-94-J-4384 


Accepted for 
publication in 
IEEE Trans. UFFC 


ESTIMATION AND CORRECTION OF ULTRASONIC 
WAVEFRONT DISTORTION USING PULSE-ECHO DATA 
RECEIVED IN A TWO-DIMENSIONAL APERTURE 


D.-L. Donald Liu* 


Robert C. Waag 


Department of Electrical Engineering 
University of Rochester 
Rochester, New York 14627 


* Currently with Siemens Medical Systems Ultrasound Group, Issaquah, WA 98029 



i 


ft 


ABSTRACT 

Pulse-echo measurements from random scattering and from a point target have been 
used to quantify transmitter beam size effects and isoplanatic patch size as well as to 
evaluate the performance of different aberration compensation techniques. Measurements 
were made using a single-element transmitter with a diameter of 1/2", 1", or 2", each fo¬ 
cused at 3". A tissue-mimicking scattering phantom or a point target was used to produce 
echoes that were received in a two-dimensional aperture synthesized by scanning a linear 
array. A specimen of abdominal wall was placed in the reception path to produce aber¬ 
ration. B-scan images were formed with no compensation, with time-shift compensation 
in the receiving aperture, and with backpropagation followed by time-shift compensation. 
The isoplanatic patch size was estimated by compensating the focus of a test point target 
with the parameters estimated for an original point target position, and observing the 
deterioration of compensation effects with increasing distance between the test and the 
original point targets. The results of the measurements using different transmitter di¬ 
ameters quantify the improvement of time-delay estimation with the increase in wavefront 
coherence that accompanies decreased transmitter beam size. For 7 specimens, the average 
isoplanatic patch size determined from a 10% increase in the —10 dB effective diameter 
was 16.7 mm in the azimuthal direction and 39.0 mm in the range direction. These sizes 
increased after backpropagation to 19.0 mm and 41.4 mm, respectively. For the 1/2", 1", 
and 2" diameter transmitters, the average contrast ratio improvement was 2.0 dB, 2.1 dB, 
and 2.8 dB, respectively, with time-shift compensation, and 2.3 dB, 2.7 dB, and 3.5 dB, 
respectively, with backpropagation of 20 mm followed by time-delay estimation and com¬ 
pensation. The investigation indicates that a tightly focused transmitter beam is necessary 
to create a scattered wavefront satisfactory for time-shift estimation, the isoplanatic patch 
is about twice as long in the range direction as in the azimuthal direction, and backprop¬ 
agation followed by time-shift compensation provides better compensation of distortion 
than time-shift compensation alone. 
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INTRODUCTION 


* 


The estimation and correction of ultrasonic wavefront distortion caused by propaga¬ 
tion through an inhomogeneous medium (as most tissues are) has been studied by academic 
and industrial researchers using various approaches. Some started from data measured with 
one-dimensional, clinical arrays having a large elevation. Some developed algorithms and 
systems with emphasis on real-time implementation. Some measured the distortion effects 
using a relatively controlled yet clinically relevant configuration. This latter approach has 
yielded transmission measurements using a spherical wave source and two-dimensional re¬ 
ceiving aperture with a small (about 1 mm 2 ) measurement spot size and has provided 
basic data that describe the characteristics of wavefront aberration caused by tissue. The 
next logical step is to apply this information to correct for tissue aberration in a pulse-echo 
setting using random scattering, similar to conditions commonly encountered in diagnostic 
imaging. A study that does so and aims to provide considerable insight into aberration 
correction in the pulse-echo mode is reported here. 

Although the existence of wavefront distortion has been known since almost the be¬ 
ginning of b-scan imaging [1,2], the need for compensation of this distortion has only 
recently become great with the advent of large, high-frequency arrays. Extensive investi¬ 
gation of aberration measurement and correction conducted in view of this need started 
about 10 years ago. Flax and O’Donnell [3] measured the time distortion based on the 
cross-correlation of waveforms received by neighboring elements of a linear array. Waag , et 
al. [4,5] reported measurements and analysis of wavefront distortion caused by propagation 
through abdominal wall and calf liver. Nock, et al. [6] proposed a ma ximum brightness 
method that was implemented in real-time. Focus compensation in these early studies 
concentrated on the correction of arrival time fluctuations. However, later studies [7-11] 
have shown that correction of arrival time fluctuation alone is not sufficient to compen¬ 
sate fully for ultrasonic distortion caused by propagation through an extended thickness 
of tissue. Measurements [7,8] indicate that the propagated wavefront contains not only 
time-shift aberration but also substantial energy level fluctuation and significant waveform 
changes. Focusing studies [9,10] show that time-shift compensation improves mainly the 
point resolution but that the sidelobes remain considerably higher than that of an ideal fo¬ 
cus. Based on these results and observations, an improved wavefront distortion model that 
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consists of a time-shifting screen placed some unknown distance away from the measure¬ 
ment aperture was proposed and employed with measured data [12]. The results obtained 
using this model with abdominal-wall data [12] and breast data [9] suggest that a better 
way to compensate for wavefront distortion is first to backpropagate the wavefront by an 
appropriate distance, and then to time-shift compensate the wavefront at that distance. 

Accurate estimation of the distortion is required to model and compensate the dis¬ 
tortion effectively. This estimation, which includes calculation of arrival time fluctuation 
from wavefront data, is a difficult task, since severe waveform distortion can invalidate 
the usual definition and estimation of arrival time. Additionally, for randomly distributed 
scatterers illuminated by a focused beam, the scattered signals are known to decorrelate 
with increasing spatial distance between the measurement spots. The degree of decorrela¬ 
tion is related to the illumination beam size, and is governed by the van Cittert-Zernike 
theorem [13,14] under the assumptions of incoherent scattering and a homogeneous prop¬ 
agation path. Another factor that influences the perceived distortion effect is the size 
of receiver elements. Simulations of a one-dimensional aperture using measurements in 
a two-dimensional aperture have shown [15] that accurate measurement of distortion re¬ 
quires small size elements to avoid significant spatial averaging over the element surface. 
All the foregoing points must be addressed for the appropriate measurement, estimation, 
and correction of wavefront distortion. 

The application of wavefront aberration correction in in vivo clinical diagnostic imag¬ 
ing requires the estimation of a distortion profile using randomly scattered waveforms 
produced by a beam of illumination. The distortion effect of the inhomogeneous medium 
may vary for different regions of the scattering medium because wavefronts originating 
from different positions propagate through different paths in the inhomogeneous medium 
and are distorted differently. Critical questions for the practical implementation of wave- 
front aberration correction are the following. How accurately can the distortion profile 
be estimated? For how large a spatial region is a single estimation suitable? How much 
improvement can be expected by distortion compensation? The answer to these questions 
naturally depends on the severity of the distortion and the compensation technique being 
used. In this paper, pulse-echo experiments and calculations designed to address the above 
questions are described to provide quantitative results about the following. 


4 



1. The effect of the transmitter beam size on the coherence of scattered wavefront and 
on the accuracy of time-delay estimation. 

2. The isoplanatic patch size, i.e., the size of a region for which the distortion effect can 
be corrected using the same set of parameters. 

3. The use of backpropagation in the pulse-echo mode using randomly scattered wave¬ 
forms, and the effectiveness of time-shift compensation with and without backpropa¬ 
gation. 

The organization of this paper is the following. First, data acquisition is described 
along with the considerations of experimental design. Next, issues in data processing, 
particularly in time-delay estimation, are discussed, and the formation of b-scan images 
with dynamic receive focusing is described for the particular geometry used in the mea¬ 
surements. Also, methods used for the evaluation of wavefront coherence, image contrast 
ratio, and point focus are described. Then, results from seven specimens of abdominal 
wall are presented for waveform similarity factor, rms arrival time fluctuation, isoplanatic 
patch size, and contrast ratio. The results are next discussed and finally conclusions are 
drawn. 
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I. DATA ACQUISITION 


A. Measurement Configuration 

The measurement apparatus consisted of a transmitter, a scattering phantom, a speci¬ 
men, and a receiving linear array, configured as shown in Fig. 1. The apparatus is immersed 
in distilled water with temperature maintained at 30° C using a heater and circulation. In 
this configuration, the transmit beam is unperturbed and only the receive beam is dis¬ 
torted by the aberration medium, for which specimens of human abdominal wall were 
used. The arrangement permits control of the transmit beam and compensation of the 
receive beam by different techniques using data measured in a two-dimensional aperture. 
Two scan movements used in the measurement are illustrated in Fig. 2. The first movement 
scans the linear array to form a two-dimensional receiving aperture. The second move¬ 
ment, analogous to conventional b-scan imaging, scans the transmitter and the receiving 
aperture together to form different scan lines. 

The transmit beam was produced using a single element transducer and the beam 
size was controlled by using transducers of different diameters, i.e., 1/2", 1", and 2", each 
focused at 3". The fixture was designed so that the focal point of each transducer was 
at the same spatial point. Due to the different electrical impedances of the transducers, 
a custom-made pulser with individual impedance matching networks for each transducer 
was used to maximize the energy radiated into the medium. 

A tissue-mimicking phantom with a 10 mm-diameter scatterer-free cylinder to simulate 
a blood vessel or a cyst was used as the scattering medium. To increase scattering signal 
strength and reduce attenuation in the phantom, the phantom was made of agar (40 grams 
per liter of water) and glass spheres of diameters around 75~90 ^m (4 grams per liter of 
agar solution). This resulted in about 10 dB greater signal strength compared to previously 
used graphite-gel phantoms. The final signal-to-noise ratio achieved with the measurement 
system averaged 21.9 dB, 22.9 dB, and 19.0 dB for seven specimens of abdominal wall using 
the 1/2", 1", and 2" transducer, respectively. 

In place of the phantom, a point target was positioned at the focal point of the trans¬ 
mitting transducer to produce a spherical scattered wavefront with the highest possible 
spatial coherence so that time delays estimated from such a wavefront could be used as a 
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standard against which time delays estimated using random scattering from the phantom 
could be compared. 

A diagram of the data acquisition system is shown in Fig. 3. A free-running oscillator, 
synchronized with the 20 MHz clock signal of the analog-to-digital (A/D) converter, trig¬ 
gered a custom-made pulser that excited the single-element transducers through impedance 
matching networks individually tuned for each transducer. The scattered signal was re¬ 
ceived by a 128-element linear array and the elements were individually accessed using 
a multiplexer. The pitch of the array is 0.72 mm and the element size in the elevation 
direction was limited to 1.00 mm to avoid excessive spatial averaging across the element 
surface. The output from the accessed element then passed through a TGC amplifier and 
a fixed gain amplifier. The signal was next A/D converted with a precision of 12 bits at a 
rate of 20 MHz. The measurement at each element was repeated 6 times (determined by 
hardware) so that data averaging could be used to reduce electronic noise and the results 
were stored in a waveform memory. After data were acquired for all elements at each 
elevation and scan position, the data were transferred to a PC via GPIB, averaged, and 
stored on a disk temporarily before being sent (via Internet) for subsequent processing on 
an IBM SP2 computer. Parameters of the measurement configuration are tabulated in 
Table I. As indicated in Table I and shown in Fig. 2, the linear array was moved over 12 
elevation positions to form the two-dimensional receiving aperture, and the transmitter 
and linear array assembly was moved over 20 positions to obtain 20 scan lines. Because 
the scan movement was mechanical, the scan step size could be chosen independent of the 
element pitch of the linear array. The size of each dataset was (2 bytes/sample) x (776 
samples/element) x (66 elements/elevation) x (12 elevations/scan) x (20 scans), or about 
25 MBytes. 

A water-path point target echo waveform and its amplitude spectrum are presented 
in Fig. 4 to show the center frequency and bandwidth of the measurement system. Typi¬ 
cal pulse-echo wavefronts obtained using a point target and using the random scattering 
phantom are shown in Fig. 5 for a representative tissue path. 

B. Measurements for Evaluating Isoplanatic Patch Size 

The point target positioned at the center of the transmit focus provided a means for 
measuring the isoplanatic patch size, i.e., the size of the scattering region for which the 
distortion effect of the inhomogeneous medium is about the same. 
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Data for the evaluation of isoplanatic patch size in the range direction were obtained 
by moving the transmitter/point-target assembly in the vertical direction manually and 
measuring the scattered wavefront in a two-dimensional aperture for each point target 
position. See Fig. 6. The topmost point target position was about 81 mm below the 
aperture, and 5 sequential positions were measured in an increment of 7 mm, yielding 5 
sets of wavefront data. 

The sets of wavefront data were independently corrected for geometric path length 
differences, after which all the wavefronts were essentially flat having only small fluctua¬ 
tions from wavefront distortion. The geometric correction was performed by estimating an 
arrival time surface, fitting a spherical surface to the arrival time surface, and shifting the 
waveforms with the fitted delays. Next, a time delay map was estimated for the wavefront 
corresponding to the topmost point target, and this was used to time shift each of the 
5 wavefronts. In this way, the wavefront corresponding to the topmost point target was 
time-shift compensated with delays estimated from itself, which is called here self compen¬ 
sation, and the other wavefronts were time-shift compensated with delays estimated for a 
different point target position, which is called here cross compensation. In cross compen¬ 
sation, a test point target wavefront is said for brevity to be compensated with distortion 
parameters estimated for the original point target. The compensated wavefronts were then 
each focused to produce an image of the point target. 

The data for the evaluation of isoplanatic patch size in the azimuthal direction was 
similarly obtained by scanning the transmitter and point-target assembly in the horizontal 
or array direction. During the scan, the specimen was stationary while the receiving array 
translated with the transmitter. This resulted in a relative shift between the specimen and 
the two-dimensional receiving aperture. In order to maintain interrogation of the same 
tissue path, a different portion of the array was used for forming the receiving aperture. 
Therefore, the azimuthal scan step size was made to be the same as the element pitch 
in the linear array, i.e., 0.72 mm, and with each scan, the portion of array used to form 
the receiving aperture was shifted by one element in the direction opposite to the scan 
movement. The number of such scans was 20. The vertical distance from the point target 
to the aperture in these measurements was about 95 mm. In processing, the 20 measured 
wavefronts were each corrected for geometric path length difference, and a time-shift map 


8 



was estimated for the first wavefront. This map was used to time-shift compensate each 
of the 20 wavefronts. The time-shift compensated wavefronts were then each focused to 
produce an image of the point target. 
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II. DATA PROCESSING 
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A. Time-Delay Estimation 

Time-delay estimation is the central issue in data processing. If time delays are not 
estimated correctly, the full potential improvement that results from time-shift compen¬ 
sation (with or without backpropagation) cannot be achieved, and the effectiveness of 
techniques that employ time-shift compensation cannot be properly assessed. S imil arly, 
for the measurement of isoplanatic patch size, if the point target focus is not optimal af¬ 
ter compensation with time delays estimated for that position, then, as the point target 
moves away, the deterioration of the compensation effect will be less noticeable, leading to 
an overestimation of the isoplanatic patch size. 


The problem of time-delay estimation for a two-dimensional array of M x N waveforms 
Si(t) for i — 1 ~ MN can be posed as an optimization problem in which the maximization 



is sought subject to a smoothness requirement on the map of Tj’s. The quantity J can be 
interpreted as the energy of the coherent sum of the aligned waveforms, and is proportional 
to the waveform similarity factor defined in Ref. [12] and the focusing criterion defined 
in Ref. [16]. The smoothness requirement on the Tj map comes from the observation 
that, since the time delay aberration is induced by a specimen of soft tissue with smooth 
structures, a two-dimensional map of arrival time is not likely to contain any abrupt jumps 
or local points that are very different from their neighboring points. This smoothness 
requirement is experimentally supported by the smoothness of point target wavefronts 
measured through tissue specimens. Although direct optimization of J is computationally 
intensive and difficult because numerous local minima exist, this general point of view was 
used to guide the development of specific time delay estimation algorithms for the study 
reported here. 


In the current study, time delays were estimated from two types of wavefront: point 
target wavefronts and random scattering wavefronts. These two types of wavefront have 
different spatial coherence properties in addition to different temporal waveform shapes. 
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The point target wavefront consists of a single pulse and is coherent throughout the aper¬ 
ture (that is, in the absence of propagation aberration, all the pulses are the same), whereas 
the scattering wavefront consists of random waveforms that lasts a long time and has a 
spatial coherence that decreases with increasing separation between points in the aperture. 
Previously, the authors have proposed two different algorithms for time-delay estimation 
in a two-dimensional aperture. These are the least-mean-square-error method [10] and 
the reference waveform method [12]. Both methods employed cross-correlation followed 
by spline interpolation to detect the peak as the basic means to find the relative time 
delay between two given waveforms. However, the reference waveform method assumes 
all the waveforms in the aperture to be the same, and so is only suitable for point target 
wavefronts. Another method known in the literature for time-delay estimation in a two- 
dimensional aperture is the row-sum method [17], but this method in its presented form 
also assumes the waveforms to be similar throughout each row and is considered to be 
suitable only for point target wavefronts. 

The least-mean-square-error method uses relative time delays for neighboring wave¬ 
forms and requires the waveforms be only locally coherent. This method is suited for 
estimating the arrival time surface for a wavefront that arises from random scattering. 
The difficulty with this method is the need to determine an arrival time when neighboring 
waveforms are poorly correlated. In the current study, the least-mean-square-error method 
was improved by using reference points to guide the search for peak positions in the cross¬ 
correlation function and using weights in the solution of the least-squares problem. The 
improved version was used as the method for time-delay estimation. 

The least-mean-square-error method for time delay estimation in a two-dimensional 
aperture has been described elsewhere [10], so only an outline is provided here. In this 
method, for a rectangular array of waveforms, relative time delays are found for neighboring 
waveforms in the horizontal, vertical, and the two diagonal directions. These relative time 
delays are then used to write a set of linear equations of the form r* — Tj = dij, where 
Ti represents the unknown time delay for waveform i and represents the given relative 
time delay between waveforms i and j. For an array of MxN waveforms, the number of 
such equations is 4 MN — 3 (M + N ) + 2, and the number of unknowns is MN — 1. (One 
waveform is used as a time reference.) Therefore, the number of equations is approximately 


11 



t 


four times the number of unknowns. These equations are solved .in the least-squares sense 

\ 

by minimizing 
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where the summations for i and j are over all the relative time delays d ^’s that have been 
given, and Wij is the weight for d^. The solution of this least-squares problem can exploit 
the simple structure of the equations and is computationally inexpensive compared to the 
calculation of all the cross-correlation functions that are needed to find the relative time 
delays. 

The estimation of relative time delay between two waveforms was based on the position 
of a peak in the cross-correlation function of the two waveforms. However, instead of 
seeking the peak position blindly, a reference point was used to guide the search for the 
peak position and to provide a weight for the obtained result. To reduce computation, 
the cross-correlation function of neighboring waveforms was calculated over only ±10 lags. 
The search in the cross-correlation function returned a peak position that was closest to 
the reference point and that position was used for the relative time delay. 

The weight of the relative time delay so obtained was calculated using the distances 
of the reference point to the two closest peaks on either side of the reference point. When 
the two distances were about equal, meaning which peak to choose was not clear, a low 
confidence was placed on the result. When one peak was much closer than the other to 
the reference point, a high confidence was placed on the result. The particular weight 
function used is log[(d 0 + e)/(<4 + e)], where d a is the larger and <4 is the smaller of the 
two distances from the reference point to the two nearest peaks, and e is a small positive 
number added to avoid singularity in the calculation. 

The reference value initially used to guide the peak search was zero for every pair, 
meaning that in the absence of any a priori knowledge, the time difference of neighboring 
waveforms was assumed to be zero. This is a reasonable initial guess so long as the 
wavefront has already been compensated for geometric path length differences. After the 
least-squares problem was solved, a new set of references was obtained using the current 
solution, i.e., Ti — Tj was used as the reference value for finding a new set of dij. With 
the newly obtained dij , the least-squares problem was solved again to yield new values for 
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Ti s. This was repeated until no dij was altered, and the final n's were used as the result 
of time-delay estimation. 


B. Dynamic Receive Focusing 

Dynamic focusing in reception was performed along the axis of the transmit beam. 
The basic idea, like that of dynamic focusing in present clinical systems, is to have the 
receive focus track the transmit wavefront. The complication in this study is that the 
transmit and receive apertures did not overlap and an angle of about 35° existed between 
the transmit axis and the normal of the receiving aperture. For tracking the transmit 
wavefront in receive focusing, the position of the transmit wavefront at each instant was 
estimated by fitting a spherical surface to the estimated arrival time surface, assuming that 
the sound speed in the scattering medium is known. A coordinate system was established 
on the surface of the two-dimensional receive aperture in which x-axis was along the 
array direction, y- axis was along the elevation direction, and 2 -axis was normal to the 
aperture surface, as diagramed in Fig. 7. The fitting algorithm used to find the center 
of a given spherical arrival time surface t(x , y) is described in an appendix. For random 
scattering, the arrival time was estimated using the portion of echoes that arose from 
the focal region of the transmission and had the highest spatial coherence. The fitted 
parameters (Xo, To, Zq) then represented the spatial origin for that portion of echoes and 
corresponded to the position of the transmit focal point. To calculate the spatial origin 
for other portions of echoes, the angle a. of the transmit axis with respect to the normal 
of the receiving aperture was used. As illustrated in Fig. 7, a distance l (At) along the 
transmission axis was determined from the focal point such that the propagation time 
along the new path is At longer than that for the path that goes through the focal point. 
The spatial origin of the wavefront that arrives At later than the focal point wavefront is 
then given by 


(x 0 (t), yo(t),z 0 (t)) = (X 0 , Yq + /(At) sina, Z 0 + /(At) cos a). 


To bring a wavefront originating from (xo(t), yo(t), zo(t)) into focus, the waveform Sj(t) 
received at (xj,j/j,0) was shifted in time by the amount 


T)i(t) =.- 

c 


\j(xo (t) - Zi) 2 + (y 0 (t) - 2 /i ) 2 + 4 (/) - \Ao(/) + 3/o(f) +*o(t) 
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This transformed the waveform s; (t) into s' (t) given by 

Si(t) = Si[t + T)i(t)]. 

Since the time-shift rji(t) is not a constant but a function of time, the waveforms were not 
only time-shifted but also stretched or compressed after dynamic focusing. To avoid unde¬ 
sirable time quantization effects, spline interpolation was used to resample the waveform 
so that rjiit) need not be quantized to be an integral number of sampling intervals. 

C. Single Transmit Image Formation 

Two kinds of images were formed with the data. The first kind is similar to the so- 
called explososcan imaging [18] in which a single transmit beam illuminates the medium 
and the received data are processed to yield an image of the medium with the given 
illumination. The second kind mimics regular b-scan imaging in which both the transmit 
and receive beams are scanned simultaneously and a line in the b-scan is obtained for each 
scan position. 

The first step common to both kinds of image formation was the estimation of a 
time-delay surface using the portion of wavefront that arises from the transmit focal point 
and has the highest spatial coherence. The signal length used for time-delay estimation is 
an important parameter. Since the signals arise from random scattering, the signals are 
themselves random. Use of a longer signal length helps overcome the effects of randomness 
on the estimated time delays. However, a longer signal extends the region in the scattering 
medium from which the wavefronts originate and the effective distortion profile of the 
specimen may consequently change. Also, the spatial coherence property of the wavefronts 
scattered from different positions may change because the transmit beam has a finite depth 
of field. Therefore, a compromise was used to determine the signal length for time-delay 
estimation. A signal length of 128 samples, corresponding to a segment length of about 
5 mm, was chosen and found to be appropriate. 

The time-delay surface was estimated by first removing an approximate geometric 
arrival time surface from the data and then applying the least-mean-square-error method 
to the shifted wavefront. The approximate shifts were then added to the least-mean-square- 
error solution to give the total time-delay surface. The reason for this two-stage process is 
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that the current implementation of the least-mean-square-error method expects a wavefront 
that does not contain significant curvature. Next, a spherical fit to the time-delay surface 
was used to find the spatial position of the focal point relative to the receiving aperture. 
The fitted spherical centers appeared to coincide with the expected transmit focal point 
positions. Then, dynamic receive focusing shifts were applied to the received wavefront 
using information about the spatial origin of the echoes obtained as described previously. 
The resulting wavefront no longer contained any spherical curvature, and was essentially 
flat in the absence of aberration. This intermediate wavefront, called here the dynamically 
shifted wavefront, plays an important role in subsequent processing. 

A single transmit image, i.e., an image formed from a single transmit beam, is s imilar 
to a snap-shot picture of the medium. In the processing, the dynamically shifted wavefront 
was windowed spatially using a Hamming window, passed through an imaging lens of fixed 
focal length, and subsequently propagated to the focal plane. The focal plane time history 
then yielded an image of the medium under the given illumination. This image can also 
be interpreted as the convolution of the transmit beam with the receive beam because 
the image is obtained by fixing the transmit beam and steering the receive beam. The 
propagation from the aperture plane to the focal plane was performed by the shift-and-add 
method using spline interpolation to implement arbitrary time-shifts of the waveforms. For 
the focal length of the imaging lens, the distance of the transmit focal point to the origin 
in the receiving aperture was used. Although the focal plane time history corresponds to a 
3-D volume, which implies that a volume image of the medium can be obtained, attention 
in the current study was limited to the time-history along the central elevation so the 
result is an image of a single cross section. 

D. B-Scan Image Formation 

B-scan images were formed without compensation (UNC), with time-shift compensa¬ 
tion (TSC), and with backpropagation followed by time-shift compensation (BP+TSC). 

To form an uncompensated b-scan image, the waveforms in the dynamically shifted 
wavefront were simply added together. This represents a focusing operation by an imaging 
lens and the result corresponds to the time history at the focal point of the imaging lens. 
The Hilbert transform [19] was used to find the analytic envelope of the waveform resulting 
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from the summation. The analytic envelopes from sequential scan lines were arranged side- 
by-side to form a b-scan image. 

♦ 

To form a time-shift compensated b-scan image, a time-delay map was estimated for 
each dynamically shifted wavefront, and the result was used to time-shift the wavefront 
at each scan position individually. After such time-shifting, the waveforms in the wave- 
front are better aligned in time, and summation results in a signal of larger amplitude. 
For wavefronts obtained with the 1/2" transmitter, time-shift compensation was also per¬ 
formed using the time-delays estimated from the corresponding 2" data. The subsequent 
processing was the same as that for the uncompensated case. 

To form an image using backpropagation followed by time-shift compensation, the 
dynamically shifted wavefront was first backpropagated by a fixed distance of 20 mm. The 
distance was fixed because of the large amount of computation needed to carry out the 
backpropagation and also because of the previous find i ng [12] that the result is relatively 
insensitive to the exact distance of backpropagation. In addition, trial computations con¬ 
firmed that the compensation effect did not change much with the present data when the 
•_/ • 

backpropagation distance was varied by a few millimeters. Since the specimen thickness 
ranged from 10 to 35 mm and a distance of 5 to 8 mm existed between the receiving aper¬ 
ture and the top surface of the specimen, backpropagation by 20 mm brought the wavefront 
about half-way into the specimen. The backpropagation was implemented using the angu¬ 
lar spectrum method [20]. In the calculation, an FFT was used to calculate the temporal 
as well as the spatial Fourier transform. Waveforms of 776 samples were zero-padded to 
1024 samples, and aperture size of 12 x 66 elements was zero-padded to 16 x 128. Since 

to backpropagate the wavefront for each scan line, the backpropagation function was com¬ 
puted for only one scan line and then stored for use with other scan lines. Evanescent wave 
components corresponding to the region / 2 +/ 2 > 1/A 2 were discarded in backpropagation 
because the propagation distance was much larger than the wavelength A. A time-delay 
map was then estimated for the backpropagated wavefront and used to time-shift the 
waveforms. The subsequent processing was the same as that for the uncompensated case. 

The flow-chart in Fig. 8 summarizes the various paths and steps taken to obtain images 
from the raw data. 


i/ a 2 - f} - 4 2 ) was needed 


the same backpropagation function of the form exp (j2irz 
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III. METHODS OF EVALUATION 

A. Evaluation of Wavefront Coherence 

Wavefront coherence was measured by the waveform similarity factor, defined for N 
waveforms Sj(t), i = 1 ~ N, by 
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in which the Tj’s are adjusted to align all the waveforms in time and to maximize the WSF. 
The value of this factor is easily shown to be always between 0 and 1, with 1 corresponding 
to the case where all waveforms are exactly alike except for a time delay and an amplitude 
scaling. The value of the WSF depends on the time delays Tj and will generally be smaller 
if the estimated t*’ s are incorrect. In the current study, the r*’s were obtained using the 
least-mean-square-error method, as summarized in Section II-A. 

B. Evaluation of Point Target Image 

Point target images were obtained using the single transmit image formation de¬ 
scribed in Section II.C. The point target was illuminated with a focused transmit beam 
and the corresponding echoes received in the two-dimensional aperture were processed un¬ 
der different conditions of compensation to produce different images. These images were 
two-dimensional, with one of the dimensions being the array (azimuthal) direction and 
the other being the time (range) direction. The evaluation of such images follows the 
procedure we used in a previous study [15]. The amplitudes of the analytic envelope were 
first projected using the maximum amplitude for each array or time position, to yield one¬ 
dimensional amplitude profiles along the time or array direction. Effective widths were 
then found from these profiles at -10, -20, and -30 dB, from the peak. The square 
root of the product of the effective widths in each direction was defined as the effective 
diameter. The effective diameter was a function of the level from the peak, and provided 
a measure of the focal spot size at different levels. Another measure was the —10 dB pe¬ 
ripheral energy ratio. This dimensionless parameter was calculated by specifying an ellipse 
whose center was at the position of maximum amplitude and whose axes were half of the 
—10 dB effective widths in each direction. The ellipse separated the image into the central 
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and the peripheral regions, and the ratio of the total energy outside the ellipse to the total 
energy inside the ellipse was the peripheral energy ratio. 

C. Evaluation of Isoplanatic Patch Size 

Using data measured as described in Section I.B, images of the test point target 
were obtained at each position after cross compensation using time-shifts derived from the 
wavefront measured with the original point target. The -10 dB effective diameter of the 
image was calculated. For self compensation, the compensation effect is the best, and the 
-10 dB effective diameter is the smallest. With increasing distance between the test and 
the original point target, the compensation effect deteriorates, and the effective diameter 
increases. When the amount of increase reached 10%, the corresponding distance between 
the original and the test point targets was taken to be half of the isoplanatic patch size. 
The reason for the factor of 2 is that the isoplanatic patch extends over both sides of the 
original point target. 

Isoplanatic patch size was evaluated for wavefronts in the aperture as well as after 
backpropagation by 20 mm. In the latter case, all the wavefronts were backpropagated 
individually first, and a set of time shifts was estimated for the backpropagated wavefront 
corresponding to the original point target. These time shifts were then applied to the other 
backpropagated wavefronts, which were then focused and evaluated as before. 

D. Evaluation of Contrast Resolution 

The evaluation of contrast resolution used b-scan images of the cyst phantom obtained 
with no compensation, with time-shift compensation in the aperture, and with backpropa¬ 
gation followed by time-shift compensation. The cyst phantom had a scatterer-free cylinder 
that was 10 mm in diameter and was positioned perpendicular to the transmit beam. To 
calculate the contrast ratio, a circular area of diameter 6 mm was manually selected in 
central part of the cyst region. Another circular area of the same diameter and 12 mm 
above the first circle was correspondingly placed in the region of uniform scattering. The 
contrast ratio was calculated as the ratio of the average energy levels in these two circles. 
The area inside the cyst region was manually selected because the alignment between the 
phantom and the scanning stages was not controlled exactly. The diameter of the circular 
area was made smaller than the known diameter of the cyst-like region in order to exclude 
relatively large amplitudes near the edge. 
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IV. RESULTS 


A total of 10 specimens were studied using the agar/glass-sphere phantom. However, 
three measurements were unusable due to misalignment of the transmit focus with the 
cyst phantom or too low SNR’s caused by specimen attenuation. Results obtained with 
the remaining seven different specimens are presented here. 

A. Waveform Similarity Factor and Time-Delay Estimation 

Representative maps of estimated arrival time fluctuation are shown in Fig. 9. These 
maps show essentially the same trend of arrival time fluctuation caused by the tissue, but 
spatial features as well as the numerical values of arrival time fluctuations estimated using 
the 2" transducer data are closer to those obtained with the point target data than those 
in the maps estimated using the 1 /2" and 1" transducer data. 

The waveform similarity for the point target and water path propagation was 0.997, 
very close to the ideal value of 1.0. The nonideality stems from the variation of the response 
characteristics of the elements in the linear array. The waveform similarity factor results 
for the seven specimens are tabulated in Table II. These results show that, at the aperture, 
the average waveform similarity factor of random scattering wavefronts was 0.261 for 1/2" 
transmitter, 0.490 for 1" transmitter, and 0.577 for 2" transmitter. The point target data 
has a much higher average waveform similarity factor at 0.935. After backpropagating the 
wavefronts by 20 mm, the average waveform similarity factor increased to 0.275, 0.513, and 
0.597, respectively, for random scattering with the 1/2", 1", and 2" transmitter, and to 
0.951 for the point target data. These results show that a point target produces a scattered 
wavefront with the highest spatial coherence, but a body of random scatterers illuminated 
with a highly focused transmit beam (such as obtained with the 2" transmitter) can be 
used to produce a wavefront of sufficient spatial coherence in the absence of isolated point 
targets. These results also show that wavefront backpropagation removes some waveform 
distortion caused by propagation through the inhomogeneous medium. 

The time-delay estimation results are summarized in Table III in which the rms value 
of the residual of fitting a spherical surface to the estimated arrival time surface is used 
to describe the fluctuation. Using the results obtained with the point target data as 
the standard, the rms arrival time fluctuation estimated using random scattering was 
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underestimated by 8.6, 5.7, and 2.9 ns for the 1/2", 1", and 2" diameter transducers, 
respectively. The corresponding average underestimation was reduced to 6.6, 3.1, and 
0.4 ns, respectively, after backpropagation by 20 mm. These results, combined with the 
above waveform similarity results, show that a tight transmit focus is needed to create a 
highly coherent scattered wavefront for the accurate measurement of time shifts, and that 
backpropagation increases the wavefront coherence and improves the accuracy of arrival 
time estimation. 

B. Point Target Images 

Point target images are shown in Fig. 10 to illustrate the compensation effects through 
a tissue path. The uncompensated image has a large spread of energy around the focal 
point. The focus is significantly improved with time-shift compensation and the sidelobe 
energy is greatly reduced. With backpropagation followed by time-shift compensation, the 
sidelobe energy is further reduced but the mainlobe remains relatively unchanged. For 
comparison, the corresponding point target image obtained through a water path is also 
shown. The image obtained with backpropagation followed by time-shift compensation 
approachs the image obtained through a water path. 

Statistics are presented in Table IV for the effective diameters of the focus at differ¬ 
ent levels from the peak as well as for the —10 dB peripheral energy ratio. The results 
show that time-shift compensation significantly improves the —10 and —20 dB effective 
diameter, but the —30 dB effective diameter after time-shift compensation (16.22 mm) is 
still much greater than the corresponding water-path value (7.79 mm). With backpropa¬ 
gation followed by time-shift compensation, the effective diameter at -30 dB is reduced 
to 8.20 mm, which is very close to the water-path value. 

C. Isoplanatic Patch Size 

A set of representative point target images used to determine the isoplanatic patch 
size is shown in Fig. 11. These images were obtained for different point target positions 
after time-shift compensation using the same time delay map estimated for the first point 
target position. The deterioration of the compensation effects as the distance between 
the test point target and the original (first) point target increases is clearly visible. For 
comparison, the point target image obtained without any compensation is also shown in 
the figure. 
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The isoplanatic patch size results are summarized in Table V. For time-shift compen¬ 
sation in the aperture, the average isoplanatic patch size was 16.7 mm in the azimuthal 
direction and 39.0 mm in the range direction. For time-shift compensation after backprop- 
agation by 20 mm, the average isoplanatic patch size increased to 19.0 mm in the azimuthal 
direction and 41.4 mm in the range direction. This indicates that backpropagation has 
effectively extended the isoplanatic patch size. In only one case (npe21) did the isoplanatic 
patch size in the range direction decrease (from 16.1 mm to 7.4 mm) with backpropagation. 
Further examination of that case reveals that the self compensated point target image has 
a —10 dB focal diameter of 2.49 mm after backpropagation, much smaller than the cor¬ 
responding value (3.29 mm) obtained with time-shift compensation alone. Therefore, the 
threshold for determining the isoplanatic patch size, which is taken as 1.1 times the —10 dB 
focal diameter with self-compensation, was correspondingly lowered. If the same threshold 
value were used to determine the isoplanatic patch size with or without backpropagation, 
the isoplanatic patch size would be 24.2 mm instead of 7.4 mm in the range direction after 
backpropagation. 

D. Images and Contrast Ratio 

A set of single transmit images are shown in Fig. 12. These images were obtained 
with the phantom illuminated by the 1/2", 1", and 2"diameter transmitters, respectively, 
and the reception is through a water path. The data used to form these images were 
from the first scan position, so the center of illumination is about 10 mm away from the 
cyst-like area. The reduction of beam size with increasing transmitter diameter is clearly 
demonstrated. In the case of the 1/2" transmitter, the illuminated region was so large 
that part of the cyst-like area appears in the lower right corner of the image. 

Another set of single transmit images are shown in Fig. 13. These images are ob¬ 
tained with the 2" transmitter through a tissue reception path (npe21) under different 
conditions of beamforming (UNC, TSC, and BP+TSC), and through a water reception 
path. These images show that the receiver beam was severely distorted by the specimen, 
was significantly improved with time-shift compensation, and was further improved with 
backpropagation followed by time-shift compensation. However, the results after compen¬ 
sation still show more energy spread compared to water path results. 
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Representative b-scan images that correspond to the single transmit images shown in 
Fig. 13 are shown in Fig. 14. The contrast ratio in this example improved from 24.0 dB 
to 26.2 dB using time-shift compensation, and to 27.5 dB using backpropagation followed 
by time-shift compensation. The water path contrast ratio was much higher at 38.0 dB. 
Statistics of the contrast ratio results for the seven tissue specimens are summarized in 
Table VI. The average improvement in contrast ratio is the smallest (2.0 dB) with time-shift 
compensation on the 1/2" transducer data, and the largest (3.5 dB) with backpropagation 
followed by time-shift compensation on the 2" transducer data. 


22 




I 


VI. DISCUSSION 

The reason for using a configuration in which the transmit beam is unperturbed is 
the need to control the transmit beam size for the purpose of this investigation. In in vivo 
b-scan imaging, however, the transmit beam as well as the receive beam will be distorted 
by the inhomogeneous medium, so any compensation of the transmit beam must be ap¬ 
plied iteratively based on current estimation of distortion. To be more specific, initially, 
a transmit focus could be formed using geometric delays, and distortion information es¬ 
timated using the scattered wavefront. Next, the distortion information can be used to 
improve the transmit focus. This would produce a tighter, transmit focus for more accu¬ 
rate estimation of distortion. Such a process could then be repeated until no improvement 
is obtained or the estimated values of distortion do not change substantially. To carry 
out this kind of iteration, however, requires a two-dimensional array with electronics that 
permits modification of time delay. If the distortion is modeled as a time-shifting screen 
placed some distance away from the aperture, then the transmit waveforms should ideally 
be individually adjustable to produce a predistorted wavefront such that subsequent prop¬ 
agation in the inhomogeneous medium will undo the distortion and a diffraction-limited 
focus will result. The predistorted wavefront would be obtained by backpropagating a 
time-shifted (but otherwise ideally focused) wavefront from the position of the hypotheti¬ 
cal phase screen to the array. When both transmit and receive beams are compensated for 
distortion, the contrast ratio improvement will be about twice as much as obtained with 
only receive beam compensation, because the effective beam pattern in b-scan imaging 
is the product of the transmit and receive beams. This implies that the contrast ratio 
improvements presented in Table VI can potentially be doubled in practice. 

The underestimation of arrival time fluctuation when the transmit beam is broad (as 
shown in Table III for the 1/2" transducer) is noteworthy. This phenomenon has been 
qualitatively observed before [3]. With a broad beam, echoes return from a large volume 
of random scatterers. Since the distortion effect for different scatterer positions is different 
(with the difference larger for scatterers farther apart), the overall distortion effect for the 
rf-signal received by a particular element is in a sense an average of the different distortion 
effects. This averaging results in the significant underestimation of arrival time fluctuation 
for a broad transmit beam. 
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A gradual variation of the perceived arrival time fluctuation with the point target 
position was observed in processing the measured wavefronts to derive the isoplanatic 
patch size. The variation consisted of changes in spatial feature details and magnitudes 
of the fluctuation. Because the relative position of the receiving array and the specimen 
remained fixed in our measurements for the isoplanatic patch size, the overall spatial 
features of the arrival time fluctuation remained unchanged. 

The size of the isoplanatic patch found in this study implies that a single distortion 
profile can be used over a reasonable region of about 20 x 40 mm 2 . However, these 
numerical values depend on many factors including the severity of distortion, the distance 
from the scattering region to the aperture, and the threshold condition used to define the 
isoplanatic patch size. Examination of Tables V reveals that the isoplanatic patch sizes 
for datasets npe21, npe24, and npe27 are relatively small. Correspondingly, the waveform 
similarity factors for these 3 sets of data are low (Table II) and the waveform arrival time 
fluctuations are high (Table III). In practice, due to the relatively slow speed of sound 
propagation in tissue, only a limited number of transmissions can be used for real-time 
image formation. Therefore, not many transmissions can be devoted to the estimation of 
the distortion profile. To circumvent this difficulty, a point of interest may be selected, 
and the distortion profile around that point may be estimated and used to compensate for 
distortion in the associated isoplanatic patch. This way the frame rate is little influenced 
and the computational requirements are also reduced. 

The expansion of the isoplanatic patch size after backpropagation merits further com¬ 
ment. If the aberration effect is indeed caused by a time-shifting screen placed at some 
distance away from the measurement aperture, then, after backpropagation by the right 
distance, the time delays required to align the waveforms will be the same no matter where 
the point target is positioned. In this case, the same set of parameters (backpropagation 
distance and arrival time fluctuations) can be used to correct for any imaging positions, 
and the isoplanatic patch size will be very large. Therefore, the increase in the isoplanatic 
patch size after backpropagation indicates that a time-shifting screen placed a certain dis¬ 
tance away is a better model for the aberration medium. This observation also indicates 
that the isoplanatic patch size is not merely a function of the aberration medium but is 
also a function of the compensation algorithm. 
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A comparison of the images in Fig. 12 shows that the beam pattern of the 2" transducer 
is only modestly narrower than that of the 1" transducer. A possible explanation is that 
the 2 "-diameter transducer has a focusing lens that is over 1 cm thicker at the edge than 
at the center. The attenuation through the lens material presumably added a spatial 
apodization to the transducer and effectively broadened the transmit beam. The small 
beamwidth difference between the 1" and 2" transducers is also reflected in the very small 
difference in contrast ratios (Table VI) achieved with these transducers. 

The results presented in Table VI show that greater improvement is achieved with 
a larger transmitting aperture. An intuitive explanation for this is that the scattered 
wavefront produced by the larger aperture is more coherent and can be aligned with greater 
accuracy. Similarly, wave backpropagation removed some of the amplitude fluctuation and 
waveform distortion developed during propagation through the specimen, and this helped 
to improve the contrast ratio further. 

The compensated images still, however, have a significantly lower contrast ratio than 
that of images obtained through a water path. For example, the average contrast ratio 
obtained with the 2" transmitter after backpropagation and time-shift compensation was 
31.3 dB, whereas the water path image in Fig. 14, also obtained using the 2" transmitter, 
had a contrast ratio of 38 dB. One reason for the incomplete compensation is that the 
backpropagation model is still an approximation of the inhomogeneous medium and does 
not completely represent all the distortion effects. Another reason is the reduced signal- 
to-noise ratio (SNR) of the received signals as a result of attenuation in the specimen. 
The average tissue-path SNR obtained with the 2" transducer was 19 dB as measured by 
digitizing the random signals from the phantom with the pulser turned on, and by digitizing 
the background noise with the pulser turned off. In a water path, however, the SNR 
obtained with the 2" transducer (npe284) was 26.6 dB. Since random noise is incoherent 
and cannot be focused, its presence increases the baseline signal energy throughout the 
image, thereby reducing the image contrast. 
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VII. CONCLUSION 


Wavefront aberration estimation and correction have been investigated using random 
scattering measured in a specially designed pulse-echo configuration in which the transmit, 
beam was unperturbed by the aberration medium. 

From the investigation, the following conclusions can be made about the accuracy of 
aberration estimation, the isoplanatic patch size, and the effectiveness of correction. 

1) A tightly focused transmit beam, such as that produced by the 2" f/1.5 transducer in 
this study, is needed to produce a sufficiently coherent scattering wavefront so that time 
delays can be accurately estimated. 

2) Accommodation of waveform shape distortion in time-delay estimation, such as in the 
modified least-mean-square-error method that incorporates a global smoothness require¬ 
ment as described in this paper, is needed to estimate the time-delay surface from a random 
scattering wavefront that has propagated through distributed inhomogeneities. 

3) The isoplanatic patch size usually decreases with increasing severity of distortion, as 
shown in this study for abdominal wall specimens. The isoplanatic patch size is also 
influenced by compensation technique, as seen in the increased isoplanatic patch size for 
backpropagation followed by time-shift compensation. 

4) Incorporation of waveform shape compensation, such as by using backpropagation prior 
to time-shift compensation as in this study, improves the contrast ratio over that obtained 
by time-shift compensation alone. 

In the absence of a priori knowledge about the propagation path, a tight transmit 
focus can only be achieved by modifying iteratively the excitation signals based on previous 
estimation of the aberration. Ideally, a two-dimensional array with programmable transmit 
waveform capability is needed for producing such a tight focus. 
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APPENDIX 


The problem is to find a spherical surface that best fits a given arrival time surface 
T{x{,yi) with an unknown offset To, 

Vi) ~ R) 4* J {Xi £o) 2 "I" (l/i Vo) 2 "H 2 2 , 

c 

where the sound speed c is assumed known, and (£ 0 , 2 / 0 , 20 ) 1S the unknown center of the ■ 
sphere. Although direct minimization of the sum of squared errors in the above approxima¬ 
tion leads to a nonlinear least-squares problem, the problem can be linearized by rewriting 
the above equation as 


c 2 [r(xi, yi) - r 0 ] 2 « (Xi - x 0 ) 2 + (2 H - 2/0 ) 2 + 4 - 


Rearranging this equation yields 


9q + X{9x 4- 2 / 1^2 + ^ 1^3 ^ bi, 


where 

n = r(xi,yi), 

90 = xl + y% + z%- (cr 0 ) 2 , 

91 = -2x 0 , 

= — 22/o, 

9 3 = 2c 2 r 0 , 

and 




Therefore, the problem can be solved by minimizing 

T: [6*0 + Xi9i + y { 9 2 + n9z - bi] 2 

i 

with respect to 9o,9\,9 2 ,9 3 . This is a linear problem, and £ 0 ,2/o, z o, and To can readily be 
found from the 9,’s. 
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Table I. Data Acquisition Parameters 


Parameter 

. Transmit Transducer diameters 
Transmit Transducer focal length 
2-D aperture in array direction 
2-D aperture in elevation direction 
Number of scan lines 
Scan line increment 
Signal center frequency 
Signal bandwidth (—6 dB) 
Sampling rate 
Sampling precision 


Value 

0.5", 1", and T 
3" 

66 x 0.72 mm 2 
12 x 1.0 mm 2 
20 

1.0 mm 
3.63 MHz 
1.42 MHz 
20 MHz 
12 bits 


Table II. Waveform Similarity Factor Statistics. The col umn under Pt are data from a 
point target, while the columns under 1/2", 1", and 2" are data from random scattering 
produced with transducers of the corresponding diameters. 


Dataset 



In Aperture 
1 / 2 " 1 " 2 


After Backpropagation 
Pt 1/2" 1" 2 


npel4 

0.977 

0.293 

0.575 

0.715 

0.964 

0.298 

0.579 

0.720 

npel5 

0.956 

0.270 

0.550 

0.674 

0.957 

0.277 

0.560 

0.685 

npel7 

0.979 

0.288 

0.571 

0.733 

0.978 

0.294 

0.575 

0.735 

npe21 

0.863 

0.224 

0.372 

0.389 

0.930 

0.251 

0.442 

0.439 

npe22 

0.935 

0.263 

0.497 

0.544 

0.957 

0.279 

0.535 

0.575 

npe24 

0.893 

0.248 

0.437 

0.474 

0.912 

0.263 

0.449 

0.490 

npe27 

0.941 

0.245 

0.429 

0.511 

0.958 

0.263 

0.449 

0.535 

avg 

0.935 

0.261 

0.490 

0.577 

0.951 

0.275 

0.513 

0.597 

sd 

0.040 

0.023 

0.073 

0.122 

0.021 

0.016 

0.059 

0.109 


Table III. Arrival Time Fluctuation Statistics. The abbreviations are the same as in Table 
II while the units are nanoseconds. 

In Aperture After Backpropagation 


Isabel 

Pt 

1 /2" 

1" 

2 " 

Pt 

1 /2" 

1" 

2 " 

npel4 

29.9 

28.6 

28.7 

28.9 

30.6 

27.5 

28.8 

29.3 

npel5 

44.3 

34.9 

40.3 

41.1 

42.4 

35.0 

41.5 

42.0 

npel7 . 

34.8 

30.4 

31.7 

33.1 

35.0 

29.4 

30.8 

32.6 

npe21 

60.9 

48.4 

52.7 

58.0 

61.8 

55.4 

60.6 

65.0 

npe22 

41.6 

30.0 

31.9 

36.0 

41.6 

33.4 

35.6 

39.0 

npe24 

62.9 

52.3 

55.5 

59.4 

61.6 

51.7 

56.5 

60.2 

npe27 

40.2 

30.6 

33.2 

39.5 

42.3 

36.5 

40.0 

44.1 

avg 

45.2 

36.6 

39.5 

42.3 

45.0 

38.4 

42.0 

44.6 

sd 

11.5 

8.9 

9.9 

11.0 

11.3 

10.0 

11.3 

12.4 
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Table IV. Point-Target Image Statistics. Listed are the average (avg) and standard devia¬ 
tion (sd) of the —10, —20, and —30 dB effective diameter in mm, and the —10 dB peripheral 
energy ratio (PER) for 7 abdominal specimens under three conditions: UNC=uncompensated, 
TSC=time-shift compensated, BP+TSC=backpropagation by 20 mm followed by time- 
shift compensation. The corresponding values for a water path are also listed. 




UNC 

TSC 

BP+TSC 

Water 

— 10 dB Dia. 

avg 

4.38 

2.40 

2.25 

2.25 


sd 

3.13 

0.08 

0.12 


-20 dB Dia. 

avg 

10.59 

4.87 

3.70 

3.38 


sd 

4.43 

2.62 

0.71 


—30 dB Dia. 

avg 

18.46 

16.22 

8.20 

7.79 


sd 

1.06 

2.42 

2.68 


-10dRPF,R aVg 

0.727 

0.297 

0.203 

0.188 


sd 

0.537 

0.123 

0.040 


a,tic Patch Size Statistics. The units 

are mm. 



In Aperture 

After Backpropagation 

-Ly CV UCIUU U A 

Azimuth 

Range 

Azimuth 

Range 

npel4 

27.4 

56.0 


27.4 

56.0 

npel5 

20.9 

51.5 


23.4 

56.0 

npel7 

27.4 

56.0 


27.4 

56.0 

npe21 

1.6 

16.1 


4.4 

7.4 

npe22 

13.0 

54.1 


15.2 

56.0 

npe24 

16.4 

15.3 


22.2 

19.6 

npe27 

9.9 

24.2 


12.9 

39.1 

avg 

16.7 

39.0 


19.0 

41.4 

sd 

8.7 

18.0 


7.9 

18.9 


Table VI. Contrast Ratio Statistics. The units are dB. Abbreviations: UNC = Uncompen¬ 
sated data, TSC = time-shift compensated data, BP+TSC = backpropagation followed by 
time-shift compensation, TSC/2" delay = time-shift compensation with delays estimated 
using the corresponding 2" transducer data. 


Dataset 

UNC 

TSC 

1/2" 

BP+ 

TSC 

TSC/ 

2" delay 

UNC 

1 " 

TSC 

BP+ 

TSC 

UNC 

2" 

TSC 

BP+ 

TSC 

npel4 

24.0 

25.0 

26.1 

24.4 

29.8 

30.6 

30.9 

30.7 

32.9 

33.0 

npel5 

21.9 

22.2 

22.7 

23.7 

28.1 

29.7 

29.6 

29.1 

30.9 

30.5 

npel7 

21.9 

23.8 

24.1 

24.9 

31.0 

31.7 

32.1 

30.5 

33.6 

34.0 

npe21 

13.7 

18.2 

18.0 

18.7 

24.2 

27.6 

28.5 

24.0 

26.2 

27.5 

npe22 

17.9 

19.7 

20.5 

20.6 

28.1 

30.8 

32.4 

28.0 

30.6 

31.9 

npe24 

16.5 

19.7 

19.5 

20.5 

26.8 

30.3 

31.2 

25.0 

29.4 

31.1 

npe27 

20.9 

21.6 

21.5 

20.9 

23.6 

26.1 

26.1 

27.5 

30.3 

31.0 

avg 

19.5 

21.5 

21.8 

22.0 

27.4 

29.5 

30.1 

27.8 

30.6 

31.3 

sd 

3.4 

2.2 

2.6 

2.2 

2.5 

1.8 

2.1 

2.4 

2.2 

1.9 
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Fig. 1. Measurement Configuration. A single-element transducer produces an unperturbed 
beam that illuminates the scattering phantom. The phantom contains a scatterer-free 
cylinder that emulates a blood vessel or a cyst-like region. Echoes are received by a 
modified linear array after propagation through a specimen of human abdominal wall. The 
transmit focal point is maintained in the same spatial position by a special transmitter 
mount (not shown) when the transmitter is changed. 





Fig. 2. Scan Movements. The linear array is scanned in the elevation direction by scan 
movement 1 to synthesize a 2-D receiving aperture. The transmitter and the Unear array 
(and, thus, the synthesized 2-D aperture) are scanned together by scan movement 2 to 
form b-scan Unes. 





Fig. 3. Diagram of Data Acquisition System. 













(xi(T2) cnpe283sa Single Pulse pk=3.59 ctr=3.63 bw=1.43 
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(xio"2) cnpe283sa Average Pulse pk=3.59 ctr=3.63 bw=1.42 
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Fig. 4. Pulse-Echo Waveforms and Corresponding Spectra obtained through a Water Path 
from a Point Target. The top row shows a waveform (left) and the corresponding spectrum 
(right) from the central array element at the central elevation. The bottom row shows the 
average pulse and its spectrum. The average pulse was obtained by adding all the pulses 
in the receiving aperture after the pulses have been aligned in time. In the waveform 
plots, the horizontal axis is sample number and the vertical axis is amplitude on a linear 
scale. The transmitter was a 2 transducer, and the point target was the rounded tip of a 
1/8 -diameter stainless steel rod positioned so that the tip echo had maximum amplitude. 

Abbreviations: pk = peak frequency, ctr = center frequency of —6 dB passband, bw = 
—6 dB bandwidth. 
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Fig. 5. Wavefronts Measured through a Tissue Path. The first three rows show wavefronts 
of scattering from the scattering phantom illuminated with the 1/2", 1", and 2" transducer, 
respectively, while the last row shows the wavefront from the point target illuminated by 
the 2 transducer. The horizontal axis is time and contains 776 samples obtained at a 
20 MHz sampling rate. The vertical axis is the array direction and spans 66 elements with 
a pitch of 0.72 mm. The rf signal amplitudes are displayed using a bipolar logarithmic gray 
scale with 40 dB dynamic range for each polarity. The dashed lines indicate the segment 
of signals used for time-delay estimation after geometric correction. 
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Fig. 6. Measurements for the Evaluation of Isoplanatic Patch Size. The scattered wavefront 
from a point target was measured for an initial (original) position (indicated by solid lines). 
The point target was then moved in the range or azimuthal directions to other (test) 
positions (indicated by dashed lines) and the scattered wavefronts were measured again. 






T(AP,0) -T(AP 0 0) = At 
P 0 =Transmit Focal Point 


Fig. 7. Estimation of the Transmit Wavefront Position along the Transmit Axis at Each 
Instant. The position of the transmit focal point Po(xo, yo» zo) was found by fitting a 
spherical surface to the arrival time surface. For echoes arriving At later than the echoes 
from the focal point, the spatial origin is l{At) away from the focal point, where £{At) 
satisfies the requirement that the propagation time along the path APiO is At longer than 
along the path APoO. The angle a between the transmit axis and the 2 -axis was assumed 
known (35° for the measurements reported here). 



Fig. 8. Flow of Data Processing and Image Formation. Raw data were first corrected 
approximately for geometric path length difference using a priori knowledge about the 
position of transmit focal point. An arrival time surface was then estimated using the 
least-mean-square-error method with echoes arising from the focal region, and a spherical 
surface was fit to the arrival time surface to determine the focal point position. Next, exact 
geometric correction and dynamic focusing were performed. A single transmit image and 
a b-scan image were then formed under three different conditions: without compensation 
(UNC), with time-shift compensation in the aperture (TSC), and with backpropagation 
followed by time-shift compensation (BP+TSC). The single transmit image used aperture 
data from a single transmit or scan line position, while the b-scan image used data from 
different transmit or scan positions. 
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Fig. 9. Tissue-Path Arrival Time Maps Estimated Using a Point Target and Using Random 
Scattering Produced by 1/2", 1", and 2" Transducers. The vertical coordinate is the array 
direction and spans 66x0.72 mm. The horizontal axis is the elevation direction and spans 
12x1.00 mm. The gray scale is linear with white corresponding to -1-100 ns and black 
corresponding to —100 ns. The standard deviation (indicated at the top of each panel) and 
the spatial features of the arrival time map estimated from random scattering approach 
those estimated using a point target as the transmit focal spot size gets smaller. The 
dataset was npe27. 
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Fig. 10. Point Target Images Obtained through a Tissue Path (Left Three Panels) and 
through a Water Path (Rightmost Panel). The tissue path images were formed without 
compensation (UNC), with time-shift compensation (TSC), and with backpropagation 
followed by time-shift compensation (BP+TSC). In each panel, the analytic envelope is 
displayed on a 40 dB logarithmic scale, and the horizontal (array) direction spans 20 mm 
while the vertical (range or time) direction consists of 776 time samples that span 29.57 mm. 
The dataset was npe22. 
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Fig. 11. Point Target Images from a Single Transmit Position for the Evaluation of Iso- 
planatic Patch Size. The panel X-Comp[0] is the image that results when the wavefront 
measured at the initial point-target position is compensated with time shifts derived from 
the same wavefront. The panels X-Comp[2], [4], [6], • • •, [18] show images that result when 
wavefronts measured with the point target shifted in the azimuthal direction by 1.44 mm, 
2.88 mm, 4.32 mm, • • •, 12.96 mm are time-shift compensated using the time shifts derived 
from X-Comp[0]. The panel Uncomp[0] shows the image that results from using the un¬ 
compensated wavefront measured at the initial point-target position. The scales are the 
same as in Fig. 10. The dataset was npe27. 
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Fig. 12. Images for a Water Path from a Single Transmit Position using Three Different 
Diameter Transducers. For the images from left to right, the corresponding transmitter 
diameter was 1/2", 1", and 2". In each panel, the analytic envelope is displayed on a 30 dB 
logarithmic scale, while the horizontal and vertical scales are the same as in Fig. 10. 
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Fig. 13. Images for a Water Path and for a Tissue-Path from a Single Transmit Position 
under Different Conditions. The 2" transducer was used as the transmitter for each image. 
The left three images were formed using tissue-path data without compensation, with 
time-shift compensation in the aperture, and with backpropagation by 20 mm followed by 
time-shift compensation, respectively. The rightmost image was formed using water-path 
data. The scales are the same as in Fig. 12. The dataset was npe21. 
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CR = 24.0 dB CR = 26.2 dB CR = 27.5 dB CR = 38.0 dB 


Fig. 14. B-Scan Images Corresponding to the Single Transmit Images in Fig. 13. The 
contrast ratio (CR) is the difference between the average energy levels for the scattering 
region (el) and for the scattering-free region (e2). The scales are the same as in Fig. 10. 
The dataset was npe21. 
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Propagation and Backpropagation 
for Ultrasonic Wavefront Design 
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Abstract —Wave backpropagation is a concept that can 
be used to calculate the excitation signals for an array with 
programmable transmit waveforms to produce a specified 
field that has no significant evanescent wave components. 
This concept can also be used to find the field at a distance 
away from an aperture based on measurements made in 
the aperture. For a uniform medium, three methods exist 
for the calculation of wave propagation and backpropaga¬ 
tion: the diffraction integral method, the angular spectrum 
method, and the shift-and-add method. The boundary con¬ 
ditions that are usually implicitly assumed by these meth¬ 
ods are analyzed, and the relationship between these meth¬ 
ods are explored. The application of the angular spectrum 
method to other kinds of boundary conditions is discussed, 
as is the relationship between wave backpropagation, phase 
conjugation, and the time-reversal mirror. Wave backprop¬ 
agation is used, as an example, to calculate the excitation 
signals for a ring transducer to produce a specified pulsatile 
plane wave with a limited spatial extent. 


I. Introduction 

T hree methods exist for calculating the acoustic 
field in a uniform medium given the vibration on a 
planar surface. They are the diffraction integral method 
[1], the angular spectrum method [2],[3], and the more em¬ 
pirical shift-and-add method [4]. All these methods have 
been widely used, yet the underlying conditions for the ap¬ 
plication of each method and the relationship between the 
methods have not been clearly addressed. For example, 
the need for the use of boundary conditions appropriate 
to the physical setting is obvious in the diffraction inte¬ 
gral method, yet the boundary conditions are not obvious 
in the angular spectrum method and in the shift-and-add 
method. This implies that some choice of the boundary 
conditions has already been made when the latter two 
methods are applied. 

An application in which the above issues arise is the 
determination of the element excitation for an array with 
programmable transmit waveforms to produce a specified 
acoustic wavefront at some position away from the array. 
This can be accomplished by the backpropagation of the 
specified wavefront to the position of each element and the 
assignment of the resulting signal as the excitation. The 
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DK45533 and HL150855, and the US Army Medical Research and 
Development Command Grant DAMD17-94-J-4384. 

The authors are with the Department of Electrical Engineer¬ 
ing, University of Rochester, Rochester, New York 14627 (e-mail: 
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procedure outlined here provides a new approach to design 
a specified wavefront. 

With a planar array, backpropagation can be performed 
straightforwardly using the angular spectrum method. 
With a curved array (like a convex or concave array), how¬ 
ever, the angular spectrum method is not efficient, and the 
way to backpropagate with the other two methods is not 
obvious. 

The purpose of this paper is to show the relationship 
between these methods as well as to answer the following 
four questions. 

1. What are the implicit boundary conditions in the 
angular spectrum method and in the shift-and-add 
method? 

2. How can backpropagation be accomplished using the 
diffraction integral method and the shift-and-add 
method? 

3. Should the signals in the shift-and-add method be dif¬ 
ferentiated with respect to time, and should a spherical 
spreading factor and an obliquity factor be included? 

4. How may calculations of propagation and backpropa¬ 
gation be implemented with curved apertures? 

The computation of the ultrasonic field produced by a 
transducer often uses the normal component of particle 
velocity [5],[6]. This approach assumes, explicitly or im¬ 
plicitly, a rigid baffle boundary condition, by using a zero 
value for the normal velocity outside of the transducer sur¬ 
face in the transducer plane. Another common approach 
is the impulse response method [7] that makes the shift- 
and-add evaluation computationally more efficient by an¬ 
alytically evaluating the diffraction integral in the time 
domain for transducers of simple geometry. This approach 
also assumes the rigid baffle boundary condition. A recent 
paper [8] confirmed the equivalence between the impulse 
response method and the angular spectrum method. How¬ 
ever, the use of a rigid baffle boundary condition may not 
be supported by experimental evidence in all cases. 

Measurements have been performed [9],[10] to investi¬ 
gate the influence of boundary conditions on the acoustic 
field produced by an ultrasonic transducer. The usual con¬ 
figuration encountered in medical ultrasonic imaging, dur¬ 
ing which the transducer array is placed on the subject’s 
skin, can be approximated by assuming the transducer to 
be placed at the surface of water. For the latter situation, 
a comparison of the water tank measurements with the¬ 
oretically computed results indicates that a zero pressure 
boundary condition or a pressure release surface (p = 0 
outside the transducer surface) is more appropriate than 
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a zero velocity boundary condition or a rigid baffle surface 
(-u — 0 outside the transducer surface) [9]. 

A rigid baffle boundary condition, however, is more fre¬ 
quently used in calculating acoustic wave propagation, in 
spite of the above noted experimental observations. The 
reasons for this may include the following: (a) the preva¬ 
lent influence of Rayleigh’s work, particularly the Rayleigh 
integral, in acoutics; (b) the applicability of a rigid baffle 
condition to some common situations, such as a speaker 
mounted on a wall; and c) the small angles usually en¬ 
countered in propagation calculations that make the use 
of the correct boundary conditions relatively unimportant. 
The current paucity of data suggests the need for further 
detailed studies, by theoretical, computational, and exper¬ 
imental means, of the fields produced by ultrasonic arrays 
in settings similar to those encountered in medical diag¬ 
nosis. Such studies are, for example, important for the de¬ 
velopment of quantitative ultrasonic imaging. 

The remainder of this paper is comprised of four sec¬ 
tions: Theory, Computations and Results, Discussion, and 
Conclusion. In the Theory, the angular spectrum method, 
the diffraction integral method, and the shift-and-acld 
method are summarized in a unified framework for cal¬ 
culation of wave propagation and backpropagation in a 
uniform medium in three-dimensional or two-dimensional 
spaces. An approximation is introduced to facilitate wave 
propagation calculations with a curved aperture. Compu¬ 
tations based on the analytic results are presented in the 
next section that includes, as an example, a demonstration 
of the use of wave backpropagation to obtain the excitation 
signals for a ring transducer to produce a spatially lim¬ 
ited plane wave. In the Discussion, the calculation of wave 
propagation for other boundary conditions is considered 
and the relationship between wave backpropagation, phase 
conjugation, and the time-reversal mirror is explored. The 
main results are summarized in the Conclusion. 


II. Theory 

Conventions for notation in this analysis are the follow¬ 
ing. Time dependence is assumed to be e~^ 1 for harmonic 
vibrations. Propagation or backpropagation is assumed to 
be along the z axis and the variable z is assumed to always 
denote a positive number. 


A . Angular Spectrum Method 

Let p(x,y,z) denote the complex temporal harmonic 
amplitude of pressure at a single frequency in a uniform 
medium. Given values of p(x, ?/, z) in the xy plane at a 
certain z, the angular spectrum is defined [2] as the Fourier 
transform of p(x, y , z) with respect to x and y. The angular 
spectrum may thus be written 

OO CO 

P(f x ,fy,z)= I f p(x,y,z)e-^f* + vMdxdy . (1) 

— OO —OO 


The relationship between P(f x -> fy> z ) an d P(fx > fy> 0 )> the 
angular spectrum at z = 0, is obtained from the fact that 
p(x,y,z) obeys the Helmholtz equation (V 2 + k 2 )p — 0, 
where k — 2tt/A is the wavenumber. The result is [2] 

P{fxi fy 5 z ) ~ P(fx> fyt fy > z ) > ( 2 ) 

where H(f x: f y ,z) is the so-called propagation function, 
which can be written explicitly as 


H(f x ,f y , z) = 

\fl + f 2 y < i/a 2 (3) 

\ e - fcs V A2 (/* 2 +4 2 )-i ; /2 + /2 > 1/A 2 . 

In the region / 2 -f / 2 < 1/A 2 , the propagation function 
has a uniform amplitude of 1, and is associated with what 
are known as homogeneous waves. In the region f 2 + f y > 
1/A 2 , the propagation function decays exponentially with 
increasing propagation distance z : and is associated with 
what are known as evanescent waves. 

In arriving at (3), a choice was made for the sign of the 
exponent. The choice was based on the physical require¬ 
ment that homogeneous waves with a time dependence of 
e~^ u;t are advanced in phase to represent a time-delay, and 
on the fact that evanescent waves decay. 

Equation (2) provides a way of calculating the field 
in an xy plane at an arbitrary z using P(f x , fy, 0) or 
p(x,y, 0). Since the calculation in (2) is a multiplication 
in the spatial-frequency domain, an inverse Fourier trans¬ 
formation to the space domain results in a convolution, 
which may be expressed 


OO OO 


p(x,y,z) — J J p(xo,yo,0)h(x 


— oo —OO 


xo,V - yo,z)dx 0 dy 0 , 



where 


h{x,y,z) 


OO OO 

J I H(f x ,fy,z)eP^ + y^df x df v (5) 


— oo —oo 


is called the impulse response function of the propagation 
process. An explicit form for h(x , y , z) can be obtained by 
examining the angular spectrum of th e free-space G reen’s 
function g(r) — e^ kr Mixr where r = \Jx 1 + y 2 + z 2 , which 
is [8],[11] 


G(f X Jy,z) = 


j A 


47T A /l-A 2 (/2 + /2) 


x e 


ifc|z| A /l-A 2 (/2+/2) 


( 6 ) 


Since apparently d G(f x , f y , z)/dz = —\H{f x ,f y) z), an 
inverse Fourier transform of this relation with respect to 


fx an d f y yields 


h(x,y,z) = -2 


d 

dz 



e jkr 

2nr 





This result is the same as that derived in [3 . 
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The backpropagation process is defined as the calcula¬ 
tion of P(fx,fy, 0) from P(f x , f y , z). From (2), the back- 
propagation function may be defined as 




e 

e 


_ e ~jA^l-A 2 (/2+/2) 
-jkzy/l-\*(f>+(f*)' f 2 + f2 < !/ A 2 

k W X2 ^+f 2 v)~\ fl + /2 > 1/A 2 . 



In the region / 2 + / 2 > 1/A 2 , H(f x ,f y , - z ) grows expo- 
nentially with z, signifying that evanescent waves get an 
exponential amplification to compensate for the exponen¬ 
tial decay they undergo in forward propagation. Because 
H(f x , fy, ~Z) is unbounded in the region / 2 + f 2 > i/a 2 , 
however, its inverse Fourier transform does not exist so 
that the impulse response function of backpropagation 
cannot be defined. 

To circumvent this problem, consider two scenarios in 
which backpropagation might be used. One is to find ex¬ 
citation signals for the production of a specified field. The 
other is to calculate the field at some distance from a mea¬ 
surement aperture. At the position of specification or mea¬ 
surement, because the spatial sampling rate is always fi¬ 
nite, the maximum spatial frequency / max that the data 
can represent is also finite, and the magnitude of angu¬ 
lar spectrum is assumed to be zero beyond / max . In these 
cases, therefore, the definition of H(f x , f y , —z) can be 
modified by letting H(f x , f y , - z ) = 0 for f x +f y > / 2 ax . 
Then, the inverse Fourier transform H(f x ,f y , —z) can be 
evaluated at least numerically. If / max > 1/A, using such 
a backpropagation function still implies the exponential 
amplification of the evanescent waves whose magnitude of 
spatial frequency is between 1/A and / max . 

Evanescent waves have negligible amplitudes if the dis¬ 
tance of propagation is more than a few wavelengths. For 
example, for an evanescent wave component with spatial 
frequency at 2/A (double the highest spatial frequency 
of nonevanescent waves), the amplitude after propagat¬ 
ing 2A will be, according to (3), reduced by e " 27rx2> <v / 3 _ 
3.5 x 1CT 10 , which is about —200 dB. This huge attenu¬ 
ation implies that to produce evanescent waves beyond a 
few wavelengths from the aperture is almost impossible, 
and to determine the existence of evanescent waves at a 
few wavelengths away is also almost impossible. In this 
case, the handling of evanescent waves becomes unimpor¬ 
tant, and an arbitrary function may be used to replace 
the growing exponential in (8) to make the inverse Fourier 
transform well defined. One particular approximation that 
leads to an analytic form of the inverse Fourier transform is 


H(f x ,fy,-z ) = H*(f x ,fy,z ) 
e -jfcz^/l-A 2 (/2 + /2) 


fx+ fy < !/A 2 

fZ + fy> 1/A 2 


( 9 ) 


which is obtained by replacing the growing exponential 
in (8) for the evanescent waves by a decaying exponential. 
The corresponding approximate backpropagation impulse 


response function is the complex conjugate of /i(x,t/, z), 


h(x,y,-z ) = -2 


d 


,-jkr 


dz V 47rr 


e -jkr 

2tt r 


- +jk 
r 


z 

r 


( 10 ) 

The results in (9) and (10) developed here for acoustic ap¬ 
plications have also been derived [12],[13] independently 
using similar reasoning and approximation for optical ap¬ 
plications. 

Since a forward propagation followed by a backpropaga¬ 
tion of equal distance should ideally restore the wavefront, 
a good way of examining h(x,y,—z) as an approximation 
of the ideal backpropagation impulse response function is 
to calculate the convolution of h(x, y, z) and h(x^y,—z) 
with respect to x and y, and to compare the result with 
the delta function, ie., to check the accuracy of 


d /e jkr \ * d f e~i kr 
dz \47rr J dz \ Att r 


« S(x)5(y) . 



B. Diffraction Integral Method 


Since p satisfies the Helmholtz equation (V 2 + k 2 )p = 0, 
the value of p at an arbitrary point r = (x,y,z) can be 
expressed as the surface integral [11] 


P(r) = 



So 


G(r| f.)&P- P (?o) SG(rTo) 


dri 0 


dn 0 


dr 0 


( 12 ) 

where So is an arbitrary surface enclosing the point r, ro 
is a point on So, and no is the outward pointing nor¬ 
mal vector of So at ro. This formula enables the calcu¬ 
lation of p(f) from simultaneous knowledge of both p(fo) 
and dp(ro)/dno on the surface So- The Green’s function 
G(r\fo) to be used in the above formula is a general solu¬ 
tion of the equation (V 2 + k 2 )G(f\ro) = — 5(f — fo), and 
consists of two terms 


'j ^ ^ 

G{r\f 0 ) = g(R) + x(r) = + x(fi) , (13) 

where R = \r — ro| is the distance between r and fo, 
g(R) is the free-space outgoing Green’s function as already 
defined, and xif) an arbitrary function satisfying the 
Helmholtz equation (V 2 + k 2 )x{r) = 0 in the volume sur¬ 
rounded by the surface So. An arbitrary choice of x(f) will 
not change the result computed by (12), but an appropri¬ 
ate choice of this function eliminates the need of simul¬ 
taneous knowledge of both p(fo) and dp(fo)/dfio on the 
surface So, thereby preventing the boundary conditions 
from being overspecified. 

Consider the case of So being comprised of the plane 
z — 0 and the hemisphere at infinity enclosing the up¬ 
per half space. The behavior of the pressure field at infin¬ 
ity is specified to satisfy the Sommerfeld radiation condi¬ 
tion [2]. Under these conditions, the formulae for comput¬ 
ing p(r) using either p(ro) alone or dp(fo)/dno alone are 
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well-known [11], 


OC OG 


P(x,y,z) = ~ J j p(x 0 ,yo,0) N 


JkR 


— oc — CO 


oo oo 


p{x,y,z ) = - 


dp(xp, yp, zp) 
dz 0 


JkR. 


dxodyo 


— oo —oo 


zo=0 


_( 15 ) 

In these expressions, R = ^/(:r — £o) 2 + (y — 2/o) 2 + • 

When pressure is specified in the z = 0 plane (known as 
the Dirichlet boundary condition), (14) applies. A special 
case of the Dirichlet boundary condition is the pressure 
release surface, in which the pressure is specified on the 
transducer surface and is assumed to be zero elsewhere in 
the transducer plane. When the normal derivative of pres¬ 
sure is specified in the z = 0 plane (known as the Neuman 
boundary condition), (15) applies. A special case of the 
Neuman boundary condition is the rigid baffle surface, in 
which the normal velocity (which is proportional to the 
normal derivative of the pressure) is specified on the trans¬ 
ducer surface and is assumed to be zero elsewhere in the 
transducer plane. A comparison of (14) with (4) and (7) 
immediately reveals that the angular spectrum computa¬ 
tion corresponds to the Dirichlet boundary condition. This 
conclusion is not surprising, since the angular spectrum 
method assumes that values of p(x, y, z) are specified in 
the plane z — 0. 

An equation that describes backpropagation and forms 
a pair with (14) is obtained using the approximate impulse 
response function of backpropagation in (10). The result is 


oo °° —jk.R 

p(x o ,y o ,0) « ^ J J p(x,y,z) N + jkj dxdy . 

— OO — OO 

(16) 

The approximation in (16) is due, as already noted, to the 
exponential decay of evanescent waves during backpropa¬ 
gation instead of exponential amplification. 

C. Propagation and Backpropagation by Shift-And-Add 

The shift-and-add method can be derived from the 
diffraction integrals by an inverse Fourier transform of (14) 
and (16) with respect to time followed by a discretization 
of the integrals in space. An approximation is customarily 
made under the assumption that the distance R is much 
larger than the wavelength A (as is usually the case). By 
omitting the term 1/R in the parentheses and denoting 
p(x, y, z) by p w (x, y, z ) to signify that the equation applies 
to a single frequency component, (14) becomes 


— oo —OO 


Pu,{ X ’V> Z ) = 


-jkz 

2n 


oo oo 


e jkR 

poj (^o? Vo ? 0) dxodyo 


—00 —00 


(IT) 


The inverse Fourier transform of (17) with respect to tem¬ 
poral frequency then yields 


dxodyo , 
(14) 


00 00 


p(x,y,z,t) = J J P( x o,yo,OA - R/c)r^ dxodyo , 


-00 —00 


. (18) 

where p(x, y , z, t) denotes the temporal pressure variation 
at point (x, y , z) and p denotes the differentiation of p with 
respect to time. The validity of the approximation leading 
to (18) requires the relation R A to hold for all the sig¬ 
nificant frequency components contained in p(x 0 , yo, 0, t). 
Subsequent discretization of (18) with respect to (x 0 ,2/o) 
leads to 


p(x,y,z,t) = —Ax 0 Ayo 

1 

X ^ ^ ^ ^ p{ x 0rm VOm 0? ^ Ft run M-sr- > (19) 


m n 


where R mn = y/(x - Xo m ) 2 + [y ~ yo-n) 2 + • Since (19) 

corresponds to time-shifting and adding signals at z = 0 
to obtain signals in the plane at z, the name shift-and-add 
is given to the method. Similarly, the spatially discretized, 
time-domain version of (16) with the assumption R » A 
can be written as 


p{ x o,y<hO,t) 


2ixc 


AxAy 


P{x p ,y q ,z,t + R pq /c )-~y > ( 20 ) 

■I f 


P q 


where R pq — 1 /(xo ~ %p) 2 + {yo ~ yq) 2 + ^• 

Both (19) and (20) come from results that were ob¬ 
tained for a planar geometry. Thus, these equations also 
only apply directly to a planar geometry, although the re¬ 
quirement is limited to the source aperture in the case 
of (19) and to the field measurement aperture in the case 
of (20). An approximate extension of these results to more 
general aperture shapes is possible, as described in the 
next section. 


D. Propagation and Backpropagation 
in Two-dimensional Space 

Assume now that the field varies only in the x and z 
directions and is constant in the y direction so the wave 
propagation is in two dimensions. The angular spectrum 
method can be applied to two-dimensional wave propaga¬ 
tion with essentially no change in form. The corresponding 
propagation function is 


H(f x ,z) = 




e jkZy/l-\ 2 fl, 

e -kz^/X 2 f 2 -l^ 


fl < i / a 2 
fl > i / a 2 


( 21 ) 
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and the corresponding approximate backpropagation func¬ 
tion is 


H(f x , ~z) = 


e -jkzy/l-X 2 fi^ f2 < 1 / X 2 

f2 > l / X 2 _ 


( 22 ) 


Using an approach similar to the one leading to (7), the 
inverse Fourier transform of H(f X: z) is found to be 


oo 


h(x,z)= I H(f x ,z)e j27ixfx df x = H^\kr) 


— OO 


(23) 


where r = a/x 2 + z 2 and 


= MZ) + jN 0 (Z) -»• y? e ' ( ^ /4) 


as £ —>• oo 
(24) 

is the zeroth-order Hankel function of the first kind. A 
check on the accuracy of the asymptotic form using the 
IMSL library [14] indicates that the relative rms error for 
the real and the imaginary parts is less than 0.1% for £ = 
hr > 100. Using the asymptotic form, the impulse response 
function of forward propagation is found to be 


h(x, z) 




^ r j(fcr- tt/4) £ 

V\r r 


(25) 


The asymptotic form of the inverse Fourier transform of 
H(f x , —z) is similarly found to be 


h(x, —z) 


VA) 


-j(kr-TT/4) 


(26) 


The above impulse response functions h(x,z) and 
h(x, — z), as well as those presented in (7) and (10) for wave 
propagation in three-dimensional space, were obtained as¬ 
suming a planar geometry. In practice, however, curved 
apertures may be encountered. To propagate from such an 
aperture or to backpropagate fields measured by such an 
aperture, an extension of these results to a curved geome¬ 
try is needed. 

A rigorous approach to treat a curved aperture is to find 
the Green’s function in the form of (13) for the particular 
geometry and then to use that Green’s function in (12). 
However, an analytic solution for x(t) i s known only for a 
few simple geometries. Therefore, this approach does not 
usually lead to uncomplicated analytic results. 

Instead, a closer examination of the planar geometry 
result using the impulse response function proves worth¬ 
while. In the following, h(x, z) in (25) for wave propagation 
in two-dimensional space will be used as an example for 
this examination. The field p(x, z) caused by a vibration 
p(x, 0) at z = 0 can be expressed using h(x, z), 


oo 


p(x, z) 



p(x o, 0)/i(# — xq, z)dx o 




f p{x o ,0)-^=e j( - kr -* /4) cos0dx o , (27) 
J v A r 


— oo 


where r = y/(x — Xq ) 2 + z 2 and cos# = z/r with # being 
the angle between the normal of the aperture and the line 
that connects the aperture point and the field point, as 
shown in Fig. 1. The factor cos# is known as the obliq¬ 
uity factor. In order to generalize this result to an arbi¬ 
trary aperture shape, the contribution of a line segment 
p(xQ,0)dx$ in the planar aperture to the field at (x,z) is 
examined. Apart from terms related to time delay, phase 
shifting, and amplitude scaling, the term cos# relates to 
the directivity of the radiation pattern of the line segment. 
Now, assume that the radiation pattern of a line segment 
in a curved aperture is little influenced by the overall shape 
of the aperture and is the same as the radiation pattern 
of a line segment in a planar aperture. (See Fig. 1.) This 
is probably a good approximation as long as the radius of 
curvature of the aperture is much larger than the wave¬ 
length (so that in the vicinity of each line segment the 
aperture is almost flat). Using this approximation, the con¬ 
tribution from individual line segments p(x o, zo)dso can be 
found using the same formula as in (27) and the total field 
is obtained by a summation. The result is 

oo 

p(x,z)& J p(x 0 ,zo)^j^e j( - kr ~ n/4) cosOdso , (28) 

— OO 


where z > z$, r = \/(x — xq ) 2 + (z — zq) 2 , dso is a line 
segment along the curved aperture, and # is the angle 
between the normal vector and the vector r. A similar 
approach was used [15] to extend the impulse response 
method (which assumes the rigid baffle boundary condi¬ 
tion) from a planar aperture to a curved aperture. Con¬ 
sider a circular aperture of radius R that spans an angle A. 
This geometry is shown in Fig. 2. The field produced by 
this aperture can be approximately computed using (28), 


p(x, z) 



p(R sin a, ~R cos a) 


1 

a/A r 


x e 


j(kr-*/4) i? 2 + r 2 - (Z 2 + Z 2 ) 


2 Rr 


Rda . (29) 


Backpropagation of a field measured along a curve can 
be computed using the same approximation with the obliq¬ 
uity factor cos # modified according to the geometry. The 
formula is 


oo 

p{x,z)&4 f p(xo,zq) e ~i( fer ~ 7r / 4 ) cos Odso , (30) 

J V A r 

— oo 

where z < z$. 

The time-domain counterparts of (28) and (30) can¬ 
not be written explicitly because the amplitude scaling 
contains a frequency-dependent factor 1/y/X = y/f /c, 
where / is temporal frequency and c is sound speed. Thus, 
the computation of a pulse wavefront propagating in two- 
dimensional space needs to be carried out for each har¬ 
monic frequency component of the pulse individually. The 
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dU, 



Plane Aperture 


cosfldUi 



Fig. 1. Transition from a planar aperture to a curved aperture. In 
a planar aperture (left), the contribution of an aperture segment 
Uodx to a field point in the normal direction is denoted by dU\. The 
contribution to a field point at the same distance but at an angle 6 
from normal is given by cosBdUi according to (27). In a gradually 
curved aperture (right), the radiation directivity of cos0 is assumed 
to be little influenced by the overall aperture shape and is used to 
compute the field contribution by a segment of the aperture. 


z 



Fig. 2. Geometry for the calculation of obliquity factor for a circu¬ 
lar aperture. The angle 9 is between the normal direction and the 
vector f that goes from the aperture point to the field point. The 
angle O' is between the z axis and the vector r. 


final result in the time domain can be obtained by adding 
the results for each harmonic frequency together. 

III. Computations and Results 

Computations were made to examine the validity of (11) 
in a particular configuration, to investigate the effects of 
the obliquity factor in three-dimensional as well as in two- 
dimensional wave propagation, and to demonstrate the use 
of backpropagation to determine element excitation of a 
ring transducer for the creation of a spatially limited plane 
wave. 

A. Examination of the Validity of (11) 

The computations were performed in the following way. 
First, complex functions h(x,y,z) and h(x,y, — z) were 
evaluated on an (x, y ) grid in the space domain. These dis¬ 
crete functions were then Fourier transformed using two- 
dimensional FFT, and the product of the transforms was 
then inverse Fourier transformed. The use of Fourier trans¬ 
form in these computations provides an efficient evaluation 
of the two-dimensional convolution in (11). 

The following parameters in the computation were cho¬ 
sen for their applicability in medical ultrasonic imaging. 
The center frequency was 3.75 MHz. The sound speed was 
1.5 mm/gs. The distance of propagation (z) was 10 mm. 
This distance was intentionally kept relatively short so 
that the impulse response function h(x, y, z) has signifi¬ 
cant values over a relatively small area and can, thus, be 
represented with fewer samples. This distance was large 
enough, however, so that evanescent waves can be safely 
ignored. An (x,y) grid spacing of 0.1 mm x 0.1 mm was 
used. This sampling interval is less than half the wave¬ 
length, which is 0.4 mm at 3.75 MHz, and satisfies the 
sampling theorem requirement. 


Test computations were performed to determine the 
necessary coverage of the xy plane, or the grid size. The ex¬ 
pected angular spectrum of h(x, y, z) has unit amplitude 

in the region ^Jff + / 2 < 1/A and decays exponentially 

beyond that. The graphical appearance of the magnitude 
of the angular spectrum is a disk. A grid size of 512 x 512 
points with the grid spacing of 0.1 mm x 0.1 mm resulted 
in four corners in the calculated angular spectrum. These 
were caused by the truncation of h(x, y, z) in space because 
only an area of 51.2 x 51.2 mm 2 was covered in this case 
while h(x,y>z) extends throughout the entire xy plane. 
When the grid size was increased to 1024 x 1024 points, 
the four corners disappeared, showing that the spatial cov¬ 
erage in this case was sufficient to make the spatial trun¬ 
cation insignificant. Zero-padding to 2048 x 2048 points in 
space was then used to avoid wraparound errors that may 
occur in the evaluation of the convolution in (11) using 
a two-dimensional FFT. (Wraparound error occurs when 
the dimensions of the expected convolution result exceed 
the dimensions of the FFT.) 

The results are illustrated in Fig. 3. These results have 
been normalized with respect to the grid spacing by mul¬ 
tiplying the FFT result and the convolution result by spa¬ 
tial sampling intervals Ax Ay. The maximum magnitude 
of h(x, y, z) in panel (a) is at x = y = 0 and equals about 
1/A z, which has the numerical value of 0.25 for A = 0.4 mm 
and z = 10 mm. In panel (b) that shows the magnitude 
of the angular spectrum of the impulse response, the di¬ 
ameter of the disk is about 1000 FFT points, whereas 
2048 FFT points corresponds to the spatial sampling fre¬ 
quency, which is 10 points per millimeter. Thus, the ra¬ 
dius is at a spatial frequency of approximately 2.5 cy/mm, 
which is the expected 1/A. The magnitude in panel c) at 
f x = f y = 0 is seen to be close to 1.0, as predicted by (3). 
Some ringing around the transition to the evanescent wave 
region is seen in the magnitude. This is a representation of 
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Fig. 3. Computations associated with the impulse response functions h(x,y,z) in (7) and h(x,y, — z) in (10). (a) Magnitude of h(x,y,z) at 
z — 10 mm. The function is sampled by 1024 x 1024 points with a spatial sampling interval of 0.1 mm and zero-padded to 2048 X 2048 points, 
(b) Magnitude of the angular spectrum of h(x,y,z). (c) Magnitude along the central horizontal line ( f y = 0) of panel (b). (d) Magnitude of 
the central 30 X 30 points that result from convolving h(x , y , z) with h(x , y, —z). (e) The real part of the convolution result along the central 
horizontal line (y = 0). 


the Gibbs phenomenon. By applying an appropriate, non- 
rectangular, window to h(x, y , z) in space domain, the ring¬ 
ing in its angular spectrum could be minimized [3], but as 
this was not a concern in the current investigation, a uni¬ 
form or rectangular window was used. The result of convo¬ 
lution of h(x, y , z) and h(x,y,—z) is presented in panel d) 
using the central 30 x 30 samples to show details of the 
result, since the whole image of 2048 x 2048 samples pro¬ 
duces important nonzero values only in the neighborhood 
of the center. The integration of the real part of the convo¬ 
lution over the 2048 x 2048 points was found to be 1.017, 
very close to the expected value of 1.0. (The imaginary 
part is essentially zero other than numerical errors, due to 
finite precision arithmetic.) The distance between the two 
minima surrounding the peak along the central horizontal 


line, shown in panel e), was found via spline interpolation 
to be about 0.66 mm. This result can be compared to an 
analytical result, obtained by taking the inverse Fourier 
transform of the disk-shaped angular spectrum (ignoring 
evanescent waves). The result is the well-known Airy func¬ 
tion of the form J\{kp)/kp with p = y/x 2 + y 2 , and the 
two minima surrounding the peak would be 1.635A apart, 
which is 0.654 mm at the wavelength of 0.4 mm. The os¬ 
cillatory tail in the plot is due to the sharp edges in the 
spectra of h(x, y, z ) and h(x, y, —z). 

A similar calculation was performed for —2 jkg(x,y,z) 
and 2 jkg*(x,y, z). These functions are good approxima- 
tions to h(x,y,z) and h(x,y, — z) with the obliquity fac¬ 
tor z/r removed, and would correspond to propagation 
and backpropagation by shift-and-add using a unit obliq- 
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Magnitudes along fy=0 



Spot Size 0.57 (mm) 



Fig. 4. Computations associated with the free-space Green’s function. 
The two panels correspond to panels c) and e) in Fig. 3. 

uity factor. The angular spectrum of -2 jkg(x,y,z) along 
f y = 0 and the result of convolution along y — 0, calcu¬ 
lated using the same parameters as used for the calcula¬ 
tions in Fig. 3, are presented in Fig. 4. As can be seen 
in Fig. 4, the magnitude of the angular spectrum, which 
equals [1-A 2 (/J+/y )] -1/2 > increases from unity to infinity 
as the spatial frequency increases from zero to 1/A. The 
singularity in the angular spectrum of g(x,y,z) is due to 
the uniform angular distribution of energy radiated from 
a point source [16]. The result of convolution is more os¬ 
cillatory in Fig. 4 than in Fig. 3, because the discontinuity 
in the angular spectrum at f% + fy = 1/A 2 is greater. The 
integration of the real part of the convolution over the 
2048 x 2048 points was found to be 1.091, farther away 
from the expected 1.0 value than the 1.017 in Fig. 3. The 
distance between the two minima surrounding the peak is 
found to be about 0.57 mm, which is 14% smaller than the 
spot size of Fig. 3. 

B. Use of Backpropagation to Determine Element 
Excitation for the Creation of a Specified Field 

The problem that initially stimulated the investigation 
presented here is to prescribe a set of waveforms that may 
be applied to a segment of the elements of a ring transducer 
to create a spatially limited plane wave at the center of the 
ring. This ring transducer was developed by a team at Uni¬ 
versity of Rochester for ultrasonic scattering and imaging 
studies. The ring transducer has a diameter of 150 mm 
and consists of 2048 elements, each able to transmit and 
receive independently with a center frequency of 2.43 MHz 
and a -6 dB bandwidth of 1.71 MHz. The transmit wave¬ 
form (pulse shape) is independently specifiable for each 
element, a unique capability not found in commercial ul¬ 
trasonic imaging systems. A diagram of the geometry is 
shown in Fig. 5. 


0.23 mm 



Fig. 5. Ring transducer and spatially limited plane wave. The ring 
is composed of 2048 elements, with the longer dimension of each el¬ 
ement being parallel to the ring axis or perpendicular to the paper 
plane. A wavefront whose amplitude a(x) along the diameter is uni¬ 
form at the central region of 50 mm is to be produced. A cosine 
roll-off extending 1 mm on each side of the central region is used to 
provide a smooth transition in amplitude. 

The field along the 150 mm diameter was discretized 
into 1000 points. This discretization meets the sampling 
requirement for waves with a frequency up to 5 MHz. The 
field was specified to have a uniform amplitude for the 
central 50 mm with an additional 1 mm cosine roll-off at 
each end, and to have zero amplitude along the remain¬ 
der of the diameter. The pulse waveform was specified to 
have a Gaussian envelope, a center frequency of 2.43 MHz, 
and a —6 dB bandwidth of 1.71 MHz. The field was thus 
represented as 

p(x,t) = a(x)e~ t2/2a r sin(27r/ 0 t) , (31) 

where a(x) is as shown in Fig. 5, /o = 2.43 MHz, and cft 
was adjusted to meet the -6 dB bandwidth requirement. 
The pulse wavefront was temporally sampled at 20 MHz. 

The basic approach was to backpropagate the specified 
field to the position of each element and prescribe the cal¬ 
culated vibration as the excitation signal for that element. 
The field produced by the prescribed excitation was then 
calculated (using forward propagation) and compared to 
the specified field. Experimental verification is beyond the 
scope of this paper (but is planned for inclusion in another 
manuscript). 

In the computations, two-dimensional wave propaga¬ 
tion was used as opposed to three-dimensional wave prop¬ 
agation for three reasons. First, putting the physical ring 
transducer aside, a check of the approximations used to ob¬ 
tain (28) and (30), ie., to see how well a backpropagation 
to a curved aperture followed by a forward propagation 
from that curved aperture restores the original wavefront, 
is of interest. Second, the far-field point of each element in 
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the elevation (height) direction, estimated using h 2 / 4A, in 
which h — 25 mm is element height and A = 0.6 mm is the 
center frequency wavelength, is 260 mm away from the 
element. This is significantly larger than the ring radius 
(75 mm) and indicates that each element is closer to a line 
source than a point source and that wave spread in the ele¬ 
vation direction is small within the range of consideration. 
Third, the ring transducer is essentially a one-dimensional 
array so that each element of the transducer can only act 
as a line transmitter or a line receiver (of finite length) 
and cannot be divided into smaller units along the eleva¬ 
tion direction. Thus, three-dimensional wave propagation 
results cannot be applied directly to the ring transducer 
system. 

The results of the two-dimensional wave propagation 
and backpropagation calculations described above are 
shown in Fig. 6. The number of elements used for the pro¬ 
duction of the spatially-limited plane wave was 768 (out 
of the total of 2048). This number was chosen based on 
tests that showed that the calculated excitation falls off 
to a negligible amplitude beyond that aperture size. To 
reduce computation, the forward propagation of the exci¬ 
tation signals was made with a coarser spatial sampling, 
in which 100 positions were employed to cover the ring 
diameter (150 mm). A small error (about 0.9% in magni¬ 
tude) was observed at the border of the spatially limited 
plane wavefront. The error is attributed to the existence 
of a nonnegligible amount of evanescent wave energy in 
the specified field. In fact, when the width of cosine roll¬ 
offs was increased from 1 mm to 5 mm, the error reduced 
dramatically to 0.03%. 

In the above computation, (29) was used to calculate 
the forward propagation from the circular aperture. For 
comparison, computations were also made to investigate 
the error that can be caused by using incorrect obliquity 
factors. Two incorrect obliquity factors were considered: 
unity and cos 0' with O' being the angle between the z axis 
and the vector r (Fig. 2). In Fig. 7, the amplitudes across 
the diameter are plotted for the originally specified field 
and for three other fields, all obtained with backpropaga¬ 
tion of the specified field to the ring followed by forward 
propagation, but using the three different obliquity factors 
during forward propagation. The error caused by the use 
of incorrect obliquity factors is clearly shown. 


IV. Discussion 

A . Nonplanar Geometry and Efficiency of Computation 

The extension of wave propagation formulae from a pla¬ 
nar aperture to curved apertures was only presented for 
two-dimensional space. However, the same kind of exten¬ 
sion applies as well to calculations of three-dimensional 
wave propagation, using either the diffraction integral 
method or the shift-and-add method, but not the angular 
spectrum method (because the last method is inherently 
based on a planar geometry). The extension is, however, 


approximate and applies only to curved apertures whose 
radius of curvature is much greater than the wavelength. 

A different approach is possible for wave propagation 
calculations associated with a nonplanar geometry. The 
approach initially propagates the field from a curved sur¬ 
face to a close-by planar surface. As long as the distance of 
propagation is small, a simple time delay or a phase shift 
(for a single frequency component) can be used to perform 
the propagation approximately. Then, the propagation or 
backpropagation starting from the planar surface can be 
performed exactly as described here. This approach is sim¬ 
ilar to the thin lens approximation used in optics. 

The numerical efficiency of the angular spectrum 
method for planar geometry deserves emphasis. For each 
single frequency, the propagation from one plane to an¬ 
other plane requires only a Fourier transform, a com¬ 
plex multiplication, and an inverse Fourier transform. 
The Fourier transforms can be implemented using a Fast 
Fourier Transform algorithm. As a result, for the cal¬ 
culation of propagation or backpropagation between two 
planes, the angular spectrum method is computationally 
the most efficient. However, this method is not applicable if 
the propagation or backpropagation starts from a curved 
boundary, in which case the other two methods may be 
employed with the above-mentioned extension. 


B. Angular Spectrum Method with 
Other Boundary Conditions 


The angular spectrum calculation in the theory sec¬ 
tion for three-dimensional wave propagation proceeds un¬ 
der the assumption that values of pressure p(x,y,z) are 
given in the aperture. The calculation requires modifica¬ 
tion when other kinds of boundary conditions are specified. 
Consider a planar surface of local reaction [11] with a uni¬ 
form impedance Z. The boundary condition imposed by 
the surface is 

1 

-p(x,y,0)+v z (x,y,0)=b(x,y) , (32) 

where v z is the normal velocity which is related to the 
normal derivative of pressure through 


v z (x,y, 0) 


1 dp(x,y,z ) 
jwp dz 



and b(x,y) is the boundary value that is usually nonzero 
on the surface of the transducer and zero elsewhere. Sub¬ 
stituting (33) into (32), taking Fourier transform with re¬ 
spect to x and y, and using the fact that the Fourier 
transform of dp/dz at z = 0 equals P(f x ffi y: 0) times 

jkyj 1 — A 2 (/2 + /y), as is evident from (2), yields 

0) = !- g(/../») _ | (34) 

where B(f X) f y ) is the Fourier transform of 6(x, y). Once 
P(fxifyi 0) is available, P(f x , f y , z) at an arbitrary z can 
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Fig. 6. Wavefront backpropagation applied to calculate ring transducer excitations for the production of a spatially limited plane wave. 
From left to right and top to bottom, respectively, are: the specified field, the calculated excitation signals, the field obtained by forward 
propagating the excitation wavefront, and the difference between the calculated and the specified field. In all the panels, horizontal axis is 
time sampled at 20 MHz. In the top-right panel, the vertical axis corresponds to element number and the display is m log scale with a 50 dB 
range for each polarity. In the other panels, the vertical axis corresponds to 100 spatial positions equally spaced along the ring diameter 
and the display is on a linear scale with gray being zero, white and black being positive and negative maximum amplitudes, respectively. 
The amplitude range is indicated in the lower right corner of each panel. 


be obtained using (2). When Z 0, the above result 
approaches that for a pressure release surface, while when 
Z —> oo, the above result approaches that for a rigid baffle 
surface. 

In the particular case of Z —> oo, B(f x ,f y ) becomes 
VzifxJy ), and (34) becomes 


P(fxJ y , 0 ) 


pCVz(fxJy) 

^1 - A2(/J + / y 2 ) 



This equation shows that a normal velocity boundary con¬ 
dition v z (x, y, 0) can be effectively transformed into a pres¬ 
sure boundary condition p(x,y,0). Multiplying (35) with 
H(f x ,fy,z) in (3), performing inverse Fourier transform 


with respect to f x and f y , and using the angular spectrum 
of free-space Green’s function in (6), yields 

p(x,y,z) = -lf 11 Vz ^ yo R ^ --dx 0 dyo , (36) 

which is the well-known Rayleigh integral in the frequency 
domain. 

C. Backpropagation, Phase Conjugation, 
and the Time-Reversal Mirror 

The relationship between wave backpropagation, phase 
conjugation, and the so-called time-reversal mirror merits 
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Fig. 7. Comparison of using correct and incorrect obliquity factors 
during forward propagation of excitation wavefront from the ring. 
From left to right and top to bottom, respectively, are: the specified 
field amplitude and the calculated field amplitudes using the correct 
obliquity factor cos# as defined in Fig. 2, using a unity obliquity 
factor, and using the obliquity factor cos#' also as defined in Fig. 2. 
The horizontal axis in these plots corresponds to 100 positions along 
the ring diameter. 


comment, since these techniques are becoming more widely 
used and the analysis of wave propagation and backprop- 
agation in this paper provides an opportunity to eluci¬ 
date the similarity and dissimilarity of these techniques. 
For the purpose of discussion, assume a real pressure field 
p{x, y, z, t) is given in an xy plane, and its temporal Fourier 
transform is denoted by y, z). Because p(x, y, z,t) is 
a real quantity, the result of phase conjugation y* (x, y, z) 
equals the Fourier transform of y(x,y,z, —t). Therefore, 
aside from differences in implementation [17], phase con¬ 
jugation and time-reversal differ only in that the former 
operates on a single frequency component whereas the lat¬ 
ter is the time-domain representation of the same opera¬ 
tion on all the frequency components. However, the time- 
domain approach has extra flexibility: the signal can be 
temporally windowed so that only a selected part of the 
signal is time-reversed. The application of phase conju¬ 
gation or the time-reversal mirror usually consists of two 
stages. First a field is measured. Then the measurement is 
either phase conjugated or time-reversed and retransmit¬ 
ted. The retransmission involves a physical propagation in 
the medium so the medium must be present. In contrast, 
wave backpropagation is a computational process in which 
no physical propagation takes place, so the medium need 
not be present. 

To examine further the relationship between backprop¬ 
agation and phase conjugation (or the time-reversal mir¬ 
ror), consider a harmonic wavefront with complex ampli¬ 
tude y(x, y, z) at a given z. The result of approximate back- 
propagation (in which evanescent waves decay instead of 
amplify with backpropagation) of p(x, y, z) to the plane z = 


0 can be computed using H* given in (9). The result is 

PBp(fxJy,0)=P(fxJy,z)H*(f x J y ,z) , (37) 

where the subscript BP denotes backpropagation. On the 
other hand, the same wavefront may be phase conjugated 
and then retransmitted to produce a wave that propagates 
back to the plane at z = 0. In this case, although the phys¬ 
ical propagation is in the — z direction, the propagated ho¬ 
mogeneous waves (/ 2 + fy < 1/A 2 ) with a time dependence 
of e~ J0jt are advanced in phase to represent a time delay, 
and the amplitude of the propagated evanescent waves de¬ 
cay [18]. Therefore, the propagation function is exactly the 
same as H{f x , f y ,z) given by (3). Since the angular spec¬ 
trum of the phase-conjugated wavefront y*(x,y,z) is 

Ppc(fxJ y ,z)= f fp*(x,y,z)e-^^+y^dxdy 

= P*(-fx,~fy,Z ) , (38) 

where the subscript PC denotes phase conjugation, the 
result of propagating a phase-conjugated wavefront in the 
—z direction is 


Ppcifx, fy, 0 ) = P*(~fx, ~fy, z)H(f x , f y , z) 

= {P(-fx,-fy,z)H*(-f x ,-f y ,z) ]* 

= ^Bp( — /xj — fyi 0) * (39) 

Interpreted in the space domain, this result shows that 

p PC (x, y, 0) = Pb p (x, y, 0) , (40) 

which, in turn, shows that the time-domain result of prop¬ 
agating a phase-conjugated wavefront is the time-reversal 
of the approximately backpropagated wavefront. 

A wavefront y(x, y, 0) that starts at plane z = 0, prop¬ 
agates to the plane at z and is measured there can, there¬ 
fore, be fully recovered by backpropagation to z = 0 if ex¬ 
ponential amplification of evanescent waves is performed 
during backpropagation. The backpropagated wavefront 
can have a resolution higher than that limited by / 2 + 
fy < l/A 2 , as demonstrated in nearfield acoustic holog¬ 
raphy where measurements were used to reconstruct the 
field [19]. On the other hand, the y(x,y, z) wavefront can 
be phase conjugated and retransmitted to propagate back 
to the plane z = 0. In the latter case, the wavefront re¬ 
sulting at z = 0 is approximately the time-reversal of the 
initial wavefront. The approximation results because the 
evanescent waves have decayed twice during this process: 
once from 0 to z and once from z to 0. Because of the de¬ 
cay, the field obtained by retransmission has a wavelength- 
limited resolution. 

A matched filter model was employed in [17] to explain 
that a time-reversal mirror can be used to achieve a point 
focus by stating that the contribution from each element 
is maximized at the same instant of time at a certain field 
point. However, this statement applies to any field point 
within a distribution of random scatterers and does not ex¬ 
plicitly identify the need for an isolated point target. The 
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need for a point target also limits the focus produced by 
time-reversal to be at the position of the point target, and 
no way apparently exists to steer the transmit beam ef¬ 
fectively to other locations (which is essential for imaging 
purpose) unless the inhomogeneous medium is modelled 
such as by a phase screen placed some distance away from 
the aperture [20]. Another impression obtained from the 
matched-filter explanation is that the result of retrans¬ 
mission contains a temporal convolution of the form of 
hi{f 0 ,T -t)* hi(r 0 ,t), where hi(r 0 ,t) is the acoustic im¬ 
pulse response from field point ro to transducer element i. 
See (9) of [17]. Therefore, pulses may become longer even 
if the acousto-electric impulse response of the elements 
hf e (t) is an ideal 8(t). The analysis presented above shows 
that this does not happen in a uniform medium if evanes¬ 
cent waves are ignored. 

To summarize, backpropagation differs conceptually 
from phase conjugation/time-reversal mirror. Backprop¬ 
agation is a computational process and can reconstruct a 
field with a resolution higher than that limited by wave¬ 
length, whereas phase conjugation/time-reversal mirror is 
a physical process and the retransmitted field always has 
a resolution that is wavelength limited. When applied for 
imaging through an inhomogeneous medium, backprop¬ 
agation assumes knowledge or requires a model of the 
medium, and can be applied for beam steering and imag¬ 
ing, whereas phase conjugation or time reversal requires 
that the medium be present and that an isolated point 
target exist to form a focus only at that point. 

V. Conclusion 

The computation of wave propagation and backprop¬ 
agation in a three-dimensional or a two-dimensional uni¬ 
form medium has been considered. Three methods, he., 
the angular spectrum method, the diffraction integral 
method, and the shift-and-add method, have been exam¬ 
ined and their relationships explored. The straightforward 
application of the angular spectrum method was shown 
to correspond to Dirichlet boundary condition. The shift- 
and-add method was shown to require a temporal differ¬ 
entiation of the signal as well as amplitude scaling and 
multiplication with an obliquity factor in order for the re¬ 
sults to be the same as obtained with the other two meth¬ 
ods. By replacing the exponential amplification of evanes¬ 
cent weaves during backpropagation with an exponential 
decay, results were obtained for backpropagation using the 
diffraction integral method and the shift-and-add method 
in addition to the angular spectrum method. Wave prop¬ 
agation in two-dimensional space was shown to differ in 
several aspects from three-dimensional wave propagation. 
In particular, no corresponding formula for shift-and-add 
method exists in the two-dimensional case. The diffraction 
integral method for two-dimensional propagation has been 
extended to apertures of arbitrary shape by changing the 
obliquity factor according to the geometry. This extension 
is based on the assumption that the aperture is smooth 


and that the contribution of each element on the aperture 
surface to a field point is little influenced by the overall 
shape of the aperture. Computations have been used to 
illustrate the above results. 
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An inverse scattering method that uses eigenfunctions of the scattering operator is presented. This 
approach provides a unified framework that encompasses eigenfunction methods of focusing and 
quantitative image reconstruction in arbitrary media. Scattered acoustic fields are described using a 
compact, normal operator. The eigenfunctions of this operator are shown to correspond to the 
far-field patterns of source distributions that are directly proportional to the position-dependent 
contrast of a scattering object. Conversely, the eigenfunctions of the scattering operator specify 
incident-wave patterns that focus on these effective source distributions. These focusing properties 
are employed in a new inverse scattering method that represents unknown scattering media using 
products of numerically calculated fields of eigenfunctions. A regularized solution to the nonlinear 
inverse scattering problem is shown to result from combinations of these products, so that the 
products comprise a natural basis for efficient and accurate reconstructions of unknown 
inhomogeneities. The corresponding linearized problem is solved analytically, resulting in a simple 
formula for the low-pass-filtered scattering potential. The linear formula is analytically equivalent to 
known filtered-backpropagation formulas for Born inversion, and, at least in the case of small 
scattering objects, has advantages of computational simplicity and efficiency. A similarly efficient 
and simple formula is derived for the nonlinear problem in which the total acoustic pressure can be 
determined based on an estimate of the medium. Computational results illustrate focusing of 
eigenfunctions on discrete and distributed scattering media, quantitative imaging of inhomogeneous 
media using products of retransmitted eigenfunctions, inverse scattering in an inhomogeneous 
background medium, and reconstructions for data corrupted by noise. © 1997 Acoustical Society 
of America. [S0001 -4966(97)02308-4] 

PACS numbers: 43.20.Fn, 43.60.Pt, 43.35.Wa, 43.80.Qf [ANN] 


INTRODUCTION 

This paper presents a new inverse scattering method that 
employs the focusing properties of certain acoustic fields ob¬ 
tained by retransmitting eigenfunctions of the scattering op¬ 
erator. 

Eigensystem decomposition of the scattering operator, 
regardless of the inversion method employed, has potential 
advantages in methods of collecting and analyzing scattering 
data. Previous work in electrical impedance tomography has 
employed eigenfunction decomposition of an operator asso¬ 
ciated with the measurement process to determine optimal 
input current patterns and quantify the achievable resolution 
of imaging systems. 1,2 These optimal inputs can also be de- 
♦ termined by iteratively retransmitting input patterns that are 

proportional to the measured scattered field. This approach is 
essentially an analog implementation of the “power 

I * method” for determining the eigenvectors of matrices. 2,3 

Likewise, the techniques of optical and acoustic phase 
conjugation 4-7 and the analogous process of time reversal 8,9 


^Current affiliation: Applied Research Laboratory, The Pennsylvania State 
University, P.O.B. 30, State College, PA 16802. 


can be understood as analog methods of computing the 
eigenfunctions of an operator associated with the phase con¬ 
jugation or time-reversal process. Simple focusing by phase 
conjugation, in which received echoes are conjugated or time 
reversed and retransmitted, is equivalent to a single iteration 
of the power method. Further iterations of this procedure 
correspond to additional steps in the power method, and thus 
converge to the most significant eigenfunction of the associ¬ 
ated operator at a rate specified by the ratio of the two largest 
eigenvalues. 3 The eigenfunctions of the “time-reversal op¬ 
erator,” whether obtained by iterative time reversal or by 
numerical diagonalization, have been previously shown to 
correspond to source distributions that can focus incident 
energy on strong, pointlike scatterers. 10,11 

Eigensystem analysis has historically played a role in 
the theory of inverse scattering for radially symmetric 
objects. 12 For these objects, separation of variables naturally 
leads to a representation of the scattering operator in terms of 
trigonometric functions. Since these eigenfunctions are the 
same for any radial scatterer, the inverse scattering problem 
could be reduced to the problem of determining the unknown 
object from the eigenvalues of the scattering operator. 
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However, before the method presented here, the focus¬ 
ing properties of eigenfunctions have not been exploited for 
quantitative reconstruction of inhomogeneous media. It is 
stated in Ref. 9 that the concept of time reversal “cannot be 
directly compared to computed tomography” or to “tech¬ 
niques that generate the image of the medium through signal 
analysis.” Although the basic principles of focusing on point 
targets using the eigenfunctions of scattering operators have 
been put forth in Ref. 10, these principles have not previ¬ 
ously been shown to apply to general distributed inhomoge¬ 
neities. Furthermore, no general imaging method has hitherto 
been based on these principles. 

The current method presents a solution to the imaging 
problem by bringing together recent results in the theory of 
focusing, diffraction tomography, and inverse problems to 
synthesize a unified framework for quantitative imaging of 
inhomogeneous media. Application of the method shows that 
focusing on distributed inhomogeneities can be achieved us¬ 
ing eigenfunctions and also provides a technique for quanti¬ 
tative imaging of discrete and distributed inhomogeneities 
using focusing properties. 

This method has several advantages over current inverse 
scattering methods. First, the eigenfunction formulation pro¬ 
vides optimal bases for reconstruction of unknown media, so 
that inversions are performed with the minimum possible 
complexity. Second, the method is applicable to any scatter¬ 
ing medium for which the total acoustic pressure associated 
with an incident plane wave can be estimated. Inverse scat¬ 
tering in inhomogeneous background media as well as itera¬ 
tive nonlinear inverse scattering can therefore be directly 
implemented. Third, part of the computation necessary for 
the inverse scattering algorithm can be performed by analog 
means using ideas from the power method. 

The present approach also provides new understanding 
about existing methods of focusing and imaging. For simple 
scattering objects, the new method presented here reduces to 
a quantitative specification of focusing similar to that ob¬ 
tained by iterative phase conjugation or time reversal. The 
eigenfunctions of scattering operators are shown not only to 
focus on pointlike scatterers, as has been previously 
shown, 10,11 but also to concentrate incident energy in the 
vicinity of general, distributed inhomogeneities. The method 
also improves on previous approaches to focusing using 
eigenfunctions in that quantitative images of medium param¬ 
eters are obtained simultaneously with optimal incident- 
wave distributions. For the case of weakly scattering objects, 
the method reduces to a simple inversion algorithm that is 
mathematically equivalent to the filtered backpropagation 
algorithm, 13-15 but is optimally tailored to the unknown scat¬ 
tering medium. The method reduces to a comparably simple 
and efficient formula for the case of weakly nonlinear in¬ 
verse scattering. 

Analysis given in Sec. I shows that eigenfunctions of 
scattering operators are equal to the acoustic fields of effec¬ 
tive source distributions that are proportional to the com¬ 
pressibility contrast of the scattering object. An inverse scat¬ 
tering method that incorporates products of retransmitted 
fields of eigenfunctions is presented. The general method is 
then employed to derive an analytic inversion formula valid 


FIG. 1. Scattering configuration. An incident plane wave traveling in the 
direction a is scattered by an inhomogeneity and the scattered field is mea¬ 
sured in the direction 0. 


under the Bom approximation as well as a simple nonlinear 
formula valid for small multiple-scattering effects. Numeri¬ 
cal implementation of these methods is presented in Sec. II. 
Numerical results shown in Sec. Ill illustrate focusing on 
discrete and distributed inhomogeneities using a few eigen¬ 
functions. Also, quantitative inverse scattering results are 
shown both within the context of a homogeneous back¬ 
ground medium and an inhomogeneous background medium. 


I. THEORY 


A. Background 

An inverse scattering method for a medium of variable 
sound speed is derived. For simplicity of exposition, the deri¬ 
vation is given for the canonical two-dimensional scattering 
configuration sketched in Fig. 1. However, with minor modi¬ 
fications, the method is applicable to arbitrary geometries 
and dimensions. 

When the incident pressure is a plane wave of unit am¬ 
plitude propagating in the direction a , so that 
p i (x) = e lka ' x , the corresponding total acoustic pressure 
p(x,a) at the position x is given by the Lippman-Schwinger 
equation 16,17 


p(x,a) = e 


ika-x 



x~y,k)q(y)p(y,a)dy. 



where G 0 (x—y) is the Green’s function for the Helmholtz 
equation in a homogeneous medium. In an unbounded two- 
dimensional medium, G 0 (x—y) is given by the Hankel func¬ 
tion (i/4)H ( 0 l \k\x~y\). xs The angle a is defined as the 
angle corresponding to the direction unit vector a , the wave 
number k is equal to 2irffc 0 where c 0 is the wave speed of 
the background medium, and / is the temporal frequency of 
the incident wave. The integral appearing in Eq. (1), as well 
as subsequent integrals in a* and y, are understood to be 
taken over the entire plane in R 2 . The scattering potential 
q is given for a medium of variable sound speed by 



The quantity within parentheses is equal, for a medium of 
constant density, to the compressibility contrast y K , as de¬ 
fined in Ref. 17. The scattering potential is assumed to be 


716 J. Acoust. Soc. Am., Vol. 102, No. 2, Pt. 1, August 1997 


T. D. Mast et al.: Eigenfunctions of the scattering operator 716 

























J 


real-valued and to be short-range, that is, the potential q 
decreases at large distances such that 

|^(x)|^C(l + |x|) _1 “ 5 , (3) 


where |x| is the magnitude of the position vector x, for some 

S>0. 


At a measurement radius r in the far field and a mea¬ 
surement angle 6 , the scattered pressure, p s ~P~Pi * is of 
the form 


p s (r, 6 ,a) = 



A( 6 >a) + o 




where A is the far-field pattern of the scattered pressure 


A(d,a)= J e ik 0 ‘ x q(x)p(x,a)dx. (5) 

The incident pressure may be more generally taken as a 
superposition of plane waves propagating in all directions, 

pM .jnay^,a. ( 6 ) 

The far-field pattern of the corresponding scattered acoustic 
pressure is then 


Af(0) = 


A(0,a)f(a)da. 



Equation (7) defines an operator A that maps an incident- 
wave distribution /( 0 ) into the corresponding far-field scat¬ 
tered pressure Af( 6 ). The operator A is related to the usual 
scattering operator S (Ref. 19) by 

S = I— - 7 —A, (8) 


where / is the identity operator. 

The operator A is compact 19 and therefore has a count¬ 
able number of discrete eigenvalues with zero as the only 
possible cluster point. In practice, only a finite number of 
eigenvalues are distinguishable from zero. Since the poten¬ 
tial q is real-valued, the scattering operator is unitary, so that 
the eigenvalues of A lie in the complex plane on the circle 
centered at -47 ri and passing through the origin. It also 
follows that A is normal (A*A=AA*, where A* is the Her- 
mitian transpose of A), so that an orthonormal basis {/,•} for 
L 2 [0,27 t] exists consisting of eigenfunctions of A. 

Since A is a normal operator, the Hermitian transpose 
A * satisfies the relation A */; = k *// » where f t is an eigen¬ 
function of A and X* is the complex conjugate of X, . The 
eigenfunctions of A therefore also satisfy the equation 

A*Af i =\\ i \ 1 f i . (9) 

Thus the functions f t also constitute a basis of eigenfunctions 
for A *A and the corresponding eigenvalues are the squared 
magnitudes of the eigenvalues of A. The operator A*A is 
essentially a far-field analog of the “time-reversal operator” 
as defined in Ref. 10. 


B. Focusing properties 

The focusing properties of A are seen by considering the 
ratio of the scattered amplitude to the incident amplitude. 
Since A is normal, the magnitude of its largest eigenvalue is 
equal to the largest possible value of this ratio for any non¬ 
zero /: 


lMl = sup 


IIA/(0)II L 2 

\\mw L 2 



where sup(-) denotes the least upper bound and II/(-)IIl 2 
denotes the root-mean-square magnitude of a square- 
integrable function. Thus the eigenfunction associated with 
the largest eigenvalue of A specifies an incident-wave distri¬ 
bution that maximizes the energy scattered to the far field. 
Other eigenfunctions also focus energy on inhomogeneities 
with an efficiency that is quantified by the associated eigen¬ 
values. 

The focusing property of eigenfunctions of A is further 
illustrated by introducing the acoustic fields of incident-wave 
distributions specified by the eigenfunctions. One may define 
retransmitted fields of an incident-wave distribution f(a) as 


E(x)=J f{a)e ika ' % da, 
F(x)= J f(a)p(x,a)da, 


( 11 ) 


where E(x) is the retransmitted field associated with the 
incident-wave distribution in a homogeneous medium and 
F(x) is the retransmitted field in a medium containing the 
inhomogeneity q(x). 

For incident-wave patterns corresponding to eigenfunc¬ 
tions that have nonzero eigenvalues, the retransmitted fields 
of Eqs. (11) can be written using Eqs. (5) and (7) in the form 


Ei(x) 



Jo(k\x~y\)F i(y)q(y)dy, 


F ,-(x) 



(p(x,e),e ik 0 ' y )Fi(y)q(y)dy. 


The brackets in Eq. (12) denote the inner product 




u( 0 )v*( 6 )d 0 , 



while the inner product appearing in the expression for E { 
was evaluated using the identity 

W=T“f e izcose dd (14) 

known as Parseval’s integral. 20 The retransmitted fields of 
Eq. (11) are thus seen to be equivalent to a weighted convo¬ 
lution of the unknown scattering potential with inner prod¬ 
ucts of acoustic fields. 

When the scattering potential q(x) is concentrated in a 
finite number of pointlike scatterers, each very small com¬ 
pared to a wavelength, Eq. (12) reduces to an expression of 
diffraction-limited focusing on each point scatterer. That is, 
for a scattering medium defined by 


717 J. Acoust. Soc. Am., Vol. 102, No. 2, Pt. 1, August 1997 


T. D. Mast et a/.: Eigenfunctions of the scattering operator 717 



M 

<?(x) = 2 Hj8(x-x j), ( 15 ) 

the retransmitted field £,(x) is 

2 77 

E,-(x)=—2 F,(x,)./ 0 (fc|x-Xj|) / u.;, (16) 

A I j 

so that in this case, the retransmitted field £,(x) is equal to a 
weighted sum of Bessel functions, each centered at the loca¬ 
tion of one of the point scatterers. These Bessel functions 
correspond to a group of diffraction-limited main lobes, cen¬ 
tered at each scatterer position Xj, with corresponding Bessel 
sidelobes that combine coherently. Thus each retransmitted 
field E t focuses to some extent on all of the individual point 
scatterers. 

The close relationship between the retransmitted fields 
of eigenfunctions and the unknown scattering potential, as 
seen in Eq. (12), is an expression of the focusing property of 
eigenfunctions. That is, since eigenfunctions of A correspond 
to incident-wave patterns that concentrate energy within the 
support of the scattering potential, they can be said to focus 
on general distributed inhomogeneities as well as pointlike 
scatterers. This idea is illustrated numerically later in this 
paper. 


C. Inverse scattering method 

Because of the focusing properties outlined above, re¬ 
transmitted fields of eigenfunctions are a useful starting point 
for inverse scattering reconstructions. A general inverse scat¬ 
tering method incorporating these ideas is outlined below. 

The starting point for this method is an expression of the 
inverse scattering problem in terms of the operator A of 
Eq. (7) and the corresponding retransmitted fields defined in 
Eq. (11): 


(Afi ,fj)= SijXj— J F i (x)EJ(x)q(x)dx, 


? ij 1 , 2 ,... . 


(17) 


The problem can be regularized by seeking the solution that 
minimizes the weighted L 2 norm 


II q\\w= j \q{x)\ 2 W{x)dx (18) 

with VF(x) an appropriate weight. For the analysis given be¬ 
low, this weight is defined as Vk(x) = (1 + |x|)^, £>0. For the 
explicit computations given later, other choices of W(x) are 
more natural. 

A solution to the minimization problem is obtained us¬ 
ing the method of Lagrange multipliers, analogous to the 
approach used in Ref. 21 for a linearized electric impedance 
tomography problem. At a minimum, the (infinite- 
dimensional) gradient of II q\\ 2 w is a linear combination of the 
gradients of the constraints in Eq. (17). The latter can be 
calculated using the two-potential formula 16 


A qi (8,a)-A q2 (6,a)= P\(x,a) 


X(q l (x)-q 2 (x))p 2 (x,0+7r)dx i 

( 19 ) 

where A q , p x , and p 2 are the scattering operators and 

the total acoustic pressures for the inhomogeneous media 
defined by q j(x) and q 2 (x ), respectively. Equation (19) 
yields the derivative 

i .A q+€ j(0,a)-Aq(0,a) f a, 


lim- 
6 —> 0 


= J p(x,a)q(x)p(x, 6+ v)dx, 

( 20 ) 


while the infinite-dimensional gradient of \\q\\ w is found from 


lim ———————-——— = 2 f q(x)q(x)W(x)dx. (21) 


6->0 


The result follows that if the potential q M {x) solves the 
regularized inverse scattering problem [minimization of the 
weighted norm from Eq. (18) under the constraint of Eq. 
(17)], q M must be of the form 


Qm( x ) ~ ^2 2 Q, m F t {x)F^x), 
Vr(xj i m 


( 22 ) 


where F* m (x), the complex conjugate of the retransmitted 
field corresponding to an incoming condition at infinity, is 
defined as 


F*(x)= j f*(a)p(x,a+ ir)da, 


(23) 


and the coefficients Q lm are the Lagrange multipliers. If the 
above gradients are taken with respect to the real and imagi¬ 
nary parts of a complex potential, Eq. (22) as stated is also 
found to be valid when the potential q M is complex. In some 
of the simplifying approximations made below, Eq. (22) will 
yield a complex potential q M even when the data are as¬ 
sumed to come from a unitary scattering operator associated 
with the real potential q. 

By substituting Eq. (22) into Eq. (17), the inverse prob¬ 
lem is reduced to the problem of finding the coefficients 
Q lm from the nonlinear system 




,=2 2 [ f 

/ m J 


F,(x)£:;(x)F / (x)F*(x) 


W(x) 


dx <2i 


• • _ i 

f J ^ 2) • • • » 


(24) 


vhere the dependence of the fields F and F* on the scatter- 
ng potential q is implicit. 

In general, the scattering potential q(x), and therefore 
;he total pressure field p(x,a), are unknown in inverse scat¬ 
hing problems. The function p(x,a) that implicitly appears 
n Eq. (24) may therefore be replaced by the best available 
estimate for the total pressure. Equation (24) can then be 
solved for the coefficients Q Xm by standard numerical tech- 
liques for solution of linear systems. 

The number of terms N can be chosen arbitrarily; how¬ 
ever, increasing N beyond the number of nonzero eigenval- 
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ues of A is of limited benefit in reconstructions. For simple 
scattering objects, q can be represented by expansions em¬ 
ploying small values of N. For instance, for an inhomogene¬ 
ity consisting of finitely many point scatterers, N comparable 
to the number of scatterers is sufficient. 

The above method simplifies further in the case of a 
weakly scattering medium, for which the total pressure p can 
be approximated by the incident pressure. In this case, taking 
the weight W(x) = l, the coefficients Q im can be evaluated 
analytically. From Eq. (22), under the Bom approximation, 
the scattering potential takes the form 


energy discrepancy is small and so is the imaginary part of 
the potential. 

The analytic solution of Eqs. (25) and (31) is equivalent 
to the well-known filtered backpropagation formula 13 " 15 and 
has the advantage of computational simplicity, as discussed 
later in this paper. Equivalence between the two formulas is 
shown by formulating an expansion of e~ lka ' x , viewed as a 
function of a, in terms of the orthonormal basis {f m (a)}. In 
view of Eqs. (11), this expansion yields the identity 

e~ ika x =2 E* m (x)f m (a). (32) 


?*(*) = 2 2 Q, m E,(x)E* m (x). (25) 

l m 

Substituting Eq. (25) into Eq. (5) gives the equation 

A(M) = X 2 Qim\ e~ ik0 ' x E l (x)E* m (x)p(x,a)dx. 

I m J 

(26) 

Replacement of p(x,a) in Eq. (26) by the incident plane 
wave e lka ‘ x , use of Eq. (11), and integration in x over R 2 
yields 

/ r\ \ 9. C TT CTT 

A(6,a)= , 2 2 2 Qim I S(0-a-0’ + a') 

fv l rn J — 77 J — TT 

XMe')r m (a')d0'da' (27) 

for 0— a not equal to 0 or tt. The double integral in Eq. (27) 
can be evaluated using the change of variables 

jt| = cos 0 f — cos a f , jr 2 = sin 0' — sin a\ (28) 

which is one-to-one when restricted to the regions a ! <6 r 
and a f > 0 f . Evaluation of the integral yields 

(277) 2 

|sin(0-a)|A(0,a) = —— 2 2 ) 

k l rn 

+f l (a+TT)f* m (d+TT)). (29) 

Equation (29) can be solved for the coefficients Q lm using 
the fact that the eigenfunctions /,( 0) are orthonormal as well 
as the reciprocity identity 16 

A(Q+ 7r,a+ 7r)=A(a, 9). (30) 

The solution is 

Qim~ — 7 f f |sin( 6— a)\A(0,a)f^(0)f m (a)da d6. 

8 TT J J 

(31) 

Equations (25) and (31) specify a solution q B to the 
linearized inverse problem. This solution is, in general, com¬ 
plex, even when the true potential q is purely real. A physi¬ 
cal way to understand why the Bom approximation yields a 
complex scattering potential for a lossless medium is to rec¬ 
ognize that this approximation neglects multiple scattering 
and thus, the resulting output energy differs from the input 
energy. The corresponding scattering operator is then no 
longer unitary, and is only physically realizable by a poten¬ 
tial with a nonzero imaginary part. For weak scattering, the 


Substituting Eq. (31) in Eq. (25) and using Eq. (32) as well 
as its conjugate gives 

q B (x)=-[\sm{a-e)\A(d,a)e ikx A»-<*)dad6, (33) 

87 T l J J 

which is the standard filtered backpropagation formula. 
Equation (33) yields the low-pass-filtered version of the true 
potential q if multiple scattering effects are negligible. The 
correct nonlinear generalization of the linearized low-pass 
filtered solution q B is the minimal L 2 (or weighted L 2 ) so¬ 
lution q M , which is of a form specified by Eq. (22). 

The inverse scattering method developed above can also 
be used with any orthonormal set of basis functions for 
L 2 [0,27r]. For instance, reconstructions can be performed 
using eigenfunctions of A for axisymmetric objects rather 
than using the eigenfunctions associated with the measured 
A. In this case, the eigenfunctions take the form 





The retransmitted fields E m can be analytically evaluated to 
be 


E m (r,<P) = ^i m e im *J m (kr), 



and the coefficients Q im for the low-pass-filtered reconstruc¬ 
tion of q are given by 

Qlm= l f A(d,a)e~ ils e ima da d6. 

(36) 

While the retransmitted fields specified by Eq. (35) are not 
ideally matched to nonaxisymmetric scattering media, they 
can be analytically evaluated and stored for use in fast re¬ 
constructions. Since these retransmitted fields are also unaf¬ 
fected by uncertainties in scattering measurements, they are 
suitable for reconstructions from data corrupted by noise. 

Finally, use of the eigenfunction method beyond linear 
inversion is demonstrated by considering the case where the 
inhomogeneous-medium retransmitted fields F can be esti¬ 
mated from a first approximation to the scattering potential 
q. One approach in this case is to solve the full system of 
equations defined by Eq. (24); however, a more numerically 
efficient correction to the Bom approximation can be ob¬ 
tained by invoking the localized nonlinear approximation in¬ 
troduced in Ref. 22 for electromagnetic scattering. This ap¬ 
proximation follows from writing the Lippman-Schwinger 
equation [Eq. (1)] in the form 
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p(x,or) = r(x)- e lka x - j (p(y,a)~p(x,a)) 

\ 

X^(y)G 0 (x-y )dy , (37) 

/ 

where the quantity T(x), called the depolarization tensor in 
electromagnetic scattering, is defined by 

If \~ l 

T(x)= 1+ q(y)G 0 (x-y)dy\ . (38) 

\ J * 

The second term in Eq. (37) is presumed to be small because 
the singularity of the Green’s function is cancelled by the 
difference term appearing in the integrand. Thus the total 
pressure may be approximated by the formula 

p(x,a)~ T(x)e ika ’ x . (39) 

The form for the scattering potential given by Eq. (22) 
then becomes 

«*(*)“ 2 QMx)E* m (x). (40) 

VK(Xj i m 

Substituting this form into Eq. (5) and using Eq. (39) 
gives 

*«U)~2 2 Qim f e-‘ k ^W(x)- , E / (x)E* m (x) 
l m J 

xr (x) 3 e ika ' x dx. (41) 

An approximate nonlinear formula for the scattering poten¬ 
tial q can be obtained by taking fk(x)^F(x) 3 . Equation (41) 
then yields the coefficients Q bn from Eq. (31) and the result¬ 
ing solution for the scattering potential is 

4m(x)~2 2 Qim —:— • (42) 

The solution of Eq. (42) is simplified by making the 

further approximation 

fst 2 -™, ( 43 ) 

which is valid for small scattering potentials. This substitu¬ 
tion results in 

<7m(x)~ 2 2 0,„ i (2£,(x)-/ 7 ;(x))£;(x). (44) 

/ m 

This nonlinear equation for the potential q M can be approxi¬ 
mately solved by using a form of the retransmitted field 
F/(x) corresponding to the low-pass-filtered potential q B or 
to another estimate of the scattering potential. 

II. COMPUTATIONAL METHODS 

The focusing and imaging methods outlined in Sec. I 
were implemented using numerically computed scattered 
fields of inhomogeneous objects. Scattering operators were 
calculated using a method due to Kirsch and Monk, 23 in 
which an inner solution of the Helmholtz equation for a me¬ 
dium of variable sound speed is matched to an outer solution 
of integral equations that implicitly satisfy the Sommerfeld 

720 J. Acoust. Soc. Am., Vol. 102, No. 2, Pt. 1, August 1997 


radiation condition. The inner solution is obtained using a 
finite-element method, while the outer integral equations are 
solved using Nystrom’s method. 

Scattering data were calculated numerically for a num¬ 
ber of incident plane waves evenly distributed over M angles 
between 0 and 277. For each incident-wave angle, the scat¬ 
tered field was computed at M far-field receiver angles be¬ 
tween 0 and 277, so that the angular sampling rate was 
MKItt) samples per radian. The number of receiver angles 
M should be chosen such that the scattered field has no sig¬ 
nificant frequency components above the Nyquist frequency 
of Ml {A it) samples per radian. This computation yields a 
discrete representation of the scattering operator A as an 
MXM matrix, A M . 

The eigenfunctions of A and their associated eigenval¬ 
ues were estimated numerically by direct computation of the 
eigenvectors and eigenvalues of A M . Retransmitted fields of 
eigenfunctions were evaluated numerically by numerical in¬ 
tegration of Eq. (11). Images of inhomogeneous objects were 
then obtained using a straightforward numerical implemen¬ 
tation of Eqs. (25) and (31). The integrals appearing in these 
equations were evaluated using corresponding discrete sum¬ 
mations of the components of A M and its eigenvectors. For 
comparison, standard diffraction tomography inversions 
were also performed by numerical integration of Eq. (33). 

Stability of the eigenfunction imaging method was 
tested by inversion of noisy data obtained by adding numeri¬ 
cally generated Gaussian white noise to the scattering matrix 
A m . The rms amplitude of the noise was specified as a frac¬ 
tion of the rms value of A M . Thus, for instance, a signal-to- 
noise ratio of 6 dB was obtained by adding noise with an rms 
amplitude one-half the rms value of A M . 

Inversions were also performed using the basis of eigen¬ 
functions corresponding to axisymmetric inhomogeneities. 
In this case, the formula of Eq. (25) was implemented nu¬ 
merically using the trigonometric basis functions defined in 
Eq. (34), the retransmitted fields given in closed form in Eq. 
(35), and the coefficients defined in Eq. (36). 

Nonlinear eigenfunction images were obtained using the 
analytic formula of Eq. (44) with the total pressure p(x,a) 
approximated by the total pressure for a medium including a 
cylinder of specified radius and compressibility contrast. 
This computation employed an exact solution for the scatter- 

i n 

ing of a plane wave by a cylinder. 

III. NUMERICAL RESULTS 

Focusing of eigenfunctions on a distributed scattering 
object is illustrated in Fig. 2. Here, the magnitudes of the 
retransmitted fields E i(x) and E 2 (x) are shown for an inho¬ 
mogeneity consisting of a weakly scattering triangle 
( 7 *.= 0.01) approximately two wavelengths in height. The 
triangle is shown in outline together with the retransmitted 
fields. The corresponding scattering operator, calculated us¬ 
ing the finite-element/Nystrom method described above, was 
represented by a matrix of size 128X 128. The retransmitted 
fields show that the significant eigenfunctions of A specify 
incident-wave patterns that concentrate energy within the 
support of the inhomogeneity. Notable is that this focused 
energy is distributed throughout the support of the triangle. 
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Implementation of the eigenfunction method in focusing 
on pointlike scatterers is illustrated in Figs. 3 and 4. These 
figures, obtained using the linearized eigenfunction method, 
show not only diffraction-limited focusing but also quantita¬ 
tive reconstructions of the associated scattering potentials. 
Figure 3 shows images made from the scattered field of two 
pointlike scatterers at locations ( — 0.5,0) and (0, —0.2), each 
of radius 0.01 and compressibility variation -0.9. The nu¬ 
merically computed scattered field was sampled at 128 
equally spaced angles for each of 128 incident-wave angles, 
so that the operator A was represented by a 128X 128 matrix. 
The wave number used was 10, so that the scatterers were 
separated by approximately one wavelength. Since, in this 
case, two eigenvalues of A were much larger than the re¬ 
maining eigenvalues, the basic reconstruction required only 
the use of two retransmitted fields. This result illustrates that, 
for an inhomogeneity consisting of finitely many pointlike 
scatterers, the present inverse scattering method provides an 
accurate reconstruction with diffraction-limited point resolu¬ 
tion using a corresponding number of eigenfunctions. 

A stability test of the eigenfunction method is illustrated 
in Fig. 4. This image shows a reconstruction of the two 
pointlike scatterers of Fig. 3 using the same scattering data 
with added Gaussian white noise for a signal-to-noise ratio 
of 3 dB. The method of reconstruction was identical to that 
used for Fig. 3. The reconstruction shown is almost indistin¬ 
guishable from the noiseless reconstruction, indicating the 
stability of the eigenfunction imaging method. 

Linear eigenfunction images of the triangular inhomoge¬ 
neity of Fig. 2 are presented in Fig. 5. These images were 
constructed using the same scattering data as that used for 
Fig. 2. The first image, obtained using five retransmitted 
fields, shows that strong focusing is achieved using only a 
few eigenfunctions of A. The entire inhomogeneity is well- 
insonified and little incident energy is transmitted outside the 
support of the inhomogeneity. The second image, obtained 
using 15 eigenfunctions, shows that the eigenfunction 
method rapidly converges to the ideal low-pass-filtered solu¬ 
tion for the scattering potential. Notable is that the eigen¬ 
function method using 15 eigenfunctions required 69.1 s of 
CPU time on a Sun SPARCstation 10, while an analogous 
image obtained using the diffraction tomography formula of 
Eq. (33), with the integrals evaluated in an analogous man¬ 
ner, required 3014.3 s. 

Eigenfunction reconstructions of a test phantom, shown 
in Figs. 6-8, illustrate application of the eigenfunction im¬ 
aging method to a larger-scale imaging problem. The phan¬ 
tom, also represented in Fig. 1, is a cylinder of compressibil¬ 
ity contrast 0.01 and diameter of 5 mm. Internal objects 
include a water-filled (cystic) region of diameter 1 mm, a 
wire of diameter 0.1 mm and compressibility contrast 
-0.5, and an internal cylinder of diameter 1 mm and com¬ 
pressibility contrast -0.01. Scattered fields were calculated 
using the methods described above, with the operator A dis¬ 
cretized as a matrix of 256X256 points. The first image 
shown in Fig. 6, obtained using the single wave number 
k— 10 has high resolution but contains ringing (Gibbs phe¬ 
nomenon) artifacts and loss of contrast in the cystic region. 
These artifacts are removed by compounding of images ob- 
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tained using five linearly spaced wave numbers, 8^£^12, 
so that the dimensionless parameter ka varied between 20 
and 30. The five-frequency image, shown in the second panel 
of Fig. 6, also shows increased point and contrast resolution 
compared to the single-frequency image. Both images shown 
in Fig. 6 were obtained using the linearized eigenfunction 
method described above, with 64 eigenfunctions of A for 
£=8, 9, and 10, 68 eigenfunctions for k = 11, and 72 eigen¬ 
functions for k— 12. 

Reconstructions of the test phantom obtained from noisy 
data are shown in Fig. 7. Gaussian white noise was added to 
the scattering data employed in Fig. 6, so that the signal-to- 
noise ratio was 6 dB at each of the frequencies employed. 
The reconstructions employed the formula of Eq. (25) and 
coefficients obtained from Eqs. (34)-(36). The numbers of 
basis functions employed were equal to the number of eigen¬ 
functions employed in Fig. 6. These results indicate the sta¬ 
bility of the method for large objects with scattering data 
severely degraded by noise. 

Nonlinear reconstructions of the same test phantom, ob¬ 
tained using Eq. (44), are presented in Fig. 8 together with 
linear reconstructions. In the nonlinear reconstructions, the 
retransmitted fields F/(x) were estimated using pressure 
fields associated with a cylinder of diameter 5 mm and com¬ 
pressibility contrast 0.01. 17 The scattering data employed 
was identical to that used in Fig. 6(b), with five linearly 
spaced wave numbers such that 20^ka^30. The number of 
eigenfunctions employed in each image were also the same 
as those used for the images in Fig. 6. The first panel shows 
the real part of the nonlinear reconstruction, taken along the 
line y = 0, together with the real part of the analogous linear 
reconstruction from Fig. 6(b). The nonlinear reconstruction 
shows improved resolution over the linear reconstruction by 
increased height of the peak associated with the internal 
wire. The second panel shows the imaginary part of the non¬ 
linear reconstruction with the corresponding linear recon¬ 
struction from Fig. 6(b). Here, the inaccuracy of the Bom 
approximation results in a significant imaginary part for the 
linear reconstruction, while the true potential is purely real. 
The nonlinear inversion shows improved quantitative accu¬ 
racy over the linear inversion by reduction of the recon¬ 
structed imaginary part. 

IV. DISCUSSION 

Our method has shown that eigenfunctions of the scat¬ 
tering operator can be employed to focus on distributed in¬ 
homogeneities as well as pointlike scatterers. However, the 
focusing on distributed inhomogeneities occurs in a different 
manner from focusing on pointlike scatterers. That is, the 
incident energy is not maximized at a single point within the 
medium. Instead, when combined according to the present 
reconstruction method, retransmitted eigenfunctions specify 
incident-wave distributions that distribute energy throughout 
the inhomogeneous region. This type of focusing, which re¬ 
sults from the eigenfunction property of maximizing the 
scattered energy, is clearly connected to imaging of the me¬ 
dium by inverse scattering. 

The quantitative inverse scattering method presented in 
this paper can considerably simplify imaging computations. 
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FIG. 2. Focusing on a distributed inhomogeneity. Magnitudes of the retransmitted fields of the two most significant eigenfunctions are shown on a linear gray 
scale with black indicating zero and white indicating maximum amplitude. The scattering object is a uniform triangle, compressibility contrast 0.01, within the 
sketched boundaries, (a) £j(x), (b) E 2 ( x). 


The computational complexity of the current method de¬ 
pends mainly on the number of significant eigenfunctions, 
which in turn depends only on the complexity of the scatter¬ 
ing medium. Furthermore, the basis for expansion of the un¬ 
known medium is determined directly from the scattering 
data. Since the basis functions are closely related to the un¬ 
known medium, reconstructions performed using this basis 
employ a minimal amount of unnecessary information. This 
property gives the present inverse scattering method advan¬ 
tages over other methods in which a fixed basis is used to 
expand the unknown medium. 24-27 These advantages are 



FIG. 3. Eigenfunction image of two pointlike scatterers, compressibility 
contrast -0.9, separated by approximately one wavelength. The image was 
obtained using retransmitted fields of the two most significant eigenfunc¬ 
tions. 
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most apparent for inhomogeneities a small number of wave- 
lengths in size. 

The present inverse scattering method also has the ad¬ 
vantage of applicability to any medium for which the back¬ 
ground pressure field can be estimated. Use of background 
pressure estimates can greatly improve accuracy over recon¬ 
structions based on simpler approximations. For instance, 
Born inversion can yield a spurious reconstructed imaginary 
part even when the true potential is real-valued; use of an 
estimated background pressure field can greatly reduce this 
error, as seen in Fig. 8. The inverse scattering method pre- 


FIG. 4. Effect of noise on eigenfunction reconstruction. The object of Fig. 3 
was reconstructed from two eigenfunctions of synthetically noised scattering 
data with a signal-to-noise ratio of 3 dB. 
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FIG. 5. Eigenfunction images of a triangle about three wavelengths in height having compressibility contrast 0.01. (a) Real part of an inversion obtained using 
retransmitted fields of five eigenfunctions, (b) Analogous image obtained using 15 eigenfunctions. 


sented here is also extensible to any background medium for 
which a pressure field can be estimated, including moving 
fluids, layered or stratified media, and enclosed or otherwise 
bounded regions. 

The imaging, focusing, and inverse scattering methods 
presented here also intrinsically take advantage of any poten¬ 
tial increase in resolution that is associated with multiple 
scattering or other higher-order effects, as long as these ef¬ 
fects are taken into account in the estimated pressure field. 
This increase in resolution is possible because the retransmit¬ 
ted fields of eigenfunctions may have desirable qualities, 
such as higher spatial-frequency components, that are asso¬ 
ciated with the presence of an inhomogeneous background. 
Such improvements in resolution have been shown experi¬ 
mentally for time-reversal focusing in a multiply scattering 
medium. 28 


The most significant eigenfunctions of A can be deter¬ 
mined experimentally through iterative retransmission of re¬ 
ceived scattered fields in a manner similar to that performed 
in Refs. 1 and 10. This implementation of the power method 3 
allows computation of the most significant eigenfunctions of 
A by analog means, which may be preferable to numerical 
computation for very large scattering objects. These eigen¬ 
functions are useful as optimal incident-wave patterns for 
inverse scattering experiments. 

The inverse scattering method presented here includes 
the assumption that the scattering potential is purely real, 
that is, the inhomogeneous medium is assumed to have zero 
absorption. Eigensystem analysis of the scattering operator 

90 TO 

A is more complicated in the presence of absorption. ’ 
However, the methods introduced here are expected to still 



FIG. 6. Eigenfunction images of a test object having background compressibility contrast 0.01, a cystic (water-filled) region, a pointlike scatterer, and an 
internal cylinder. The images are shown on a logarithmic scale with a 40 dB dynamic range, (a) Real part of inversion obtained for a wave number such that 
ka~ 25. (b) Analogous image obtained using five wavenumbers such that 20<ka<30. 
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FIG. 7. Images of test object obtained using noised data with 6 dB signal-to-noise ratio. The images are shown on a logarithmic scale with a 40 dB dynamic 
range, (a) Real part of inversion obtained for a wave number such that ka = 25. (b) Analogous image obtained using five wave numbers such that 20<ka 

<30. 



be useful for media such as biological tissue when absorption 
is finite but small. 

A disadvantage of the inverse scattering method as cur¬ 
rently implemented is that nonlinear inversion requires an 
accurate specification of the background acoustic field. This 
disadvantage is not unique to the eigenfunction method, but 
is a common feature of most current inversion methods. Re¬ 
cent theoretically exact methods, while not limited in this 
manner, 31,32 have not yet been implemented numerically. 
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FIG. 8. Cross sections of real and imaginary parts of test object reconstruc¬ 
tions using data from five wave numbers. Linear reconstructions were ob¬ 
tained using retransmitted fields in a homogeneous medium, while nonlinear 
reconstructions were obtained using retransmitted fields in a homogeneous 
medium and in a background cylinder of compressibility contrast 0.01. (a) 
Cross sections of real parts, (b) Cross sections of imaginary parts. 
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The general, nonlinear version of the eigenfunction 
imaging approach, as defined in Eqs. (24) and (44) and illus¬ 
trated in Fig. 8, has obvious application to iterative recon¬ 
struction of unknown inhomogeneities. In such reconstruc¬ 
tions, the total pressure field can be estimated at each 
iteration from a numerical solution of the direct scattering 
problem for the currently estimated inhomogeneity, and an 
eigenfunction inversion can be performed using this pressure 
field. This procedure can then be repeated to obtain im¬ 
proved estimates of the scattering potential until convergence 
is achieved. 

V. CONCLUSION 

A method for focusing and imaging using scattered 
fields has been presented. The method outlined here makes 
use of the physical properties of scattering operators by using 
their eigenfunctions as incident-wave patterns. The eigen¬ 
functions have been shown to provide optimal focusing on 
pointlike and distributed scattering objects. 

The inverse scattering scheme presented exploits the fo¬ 
cusing properties of eigenfunctions as well as recent analytic 
results to obtain a robust and efficient means of quantita¬ 
tively reconstructing unknown scattering media. This new 
method has a complexity dependent only on the size and 
complexity of the scattering medium. Particular cases of the 
method provide improved implementations of eigenfunction 
focusing and filtered backpropagation. The method can be 
implemented for any background medium for which the total 
acoustic pressure field can be estimated. 

Another particular case of the presented method results 
in a nonlinear inverse scattering formula that yields a solu¬ 
tion for the scattering potential q in terms of retransmitted 
fields of eigenfunctions in the scattering medium and in the 
background medium. This formula has been demonstrated to 
yield improvement in accuracy and resolution over Bom in¬ 
version. 
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The ideas reported here are expected to be useful in 
further studies of inverse scattering, adaptive focusing, and 
ultrasonic imaging. The eigenfunctions of the scattering op¬ 
erator A , whether determined by iterative retransmission or 
by numerical diagonalization, may be used to focus through 
inhomogeneous media and to determine optimal incident- 
wave patterns for inverse scattering experiments. Also, the 
products of fields of eigenfunctions are expected to form a 
useful basis for expansion of unknown scattering media in 
iterative reconstruction algorithms. 
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